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1 Introduction. 

A study has been conducted by the University of South Florida 
to evaluate elements of Advanced Solar Dynamic power generation. 
This effort has been designed to determine the thermal performance 
of basic receiver concepts and the effects of various phase change 
materials. The study was funded by the NASA-Lewis Research Center 
under contract NAG3-851 and conducted during the period of January 
1988 through December 1988 at the University of South Florida, 
Tampa, Florida. 

Associated with a significant increase in space exploration, 
scientific research, commercial development and potential military 
weapons systems, comes the need for increased capacity, space 
operable, electrical power sources. Photovoltaic power generation 
is now a well developed technology successfully utilized in numerous 
commercial and military satellites. The technology is especially 
attractive for space applications in terms of mechanical simplicity, 
utilizing no moving parts and easing problems of deployment and 
maintenance. Their major limitation is relatively low energy 
efficiency, typically in the range of 8 percent. This is well below 
the levels encountered in many relatively simple terrestrial systems 
which regularly achieve efficiencies of 20 to 25 percent. 

As power requirements in space increase, the size penalties 
placed on design by the photovoltaic technology become less 
acceptable. Whereas current satellites have generally required 
less than 10 kWe, missions planned within the next 10 years [1] 
require one to two orders of magnitude increase. Advanced phases 
of the space station incorporate design plans for a number of solar 
dynamic cycle electrical generators providing a total power level 
of between 300 and 1000 kWe. Other unclassified large scale missions 
include a GEO communications platform, an air/ocean traffic control 
system, a GEO payload delivery system and a micro g materials 
processing system. With a clear need established for an advanced, 
light weight, increased capacity power generation system, NASA has 
undertaken the development of a space solar dynamic technology. 
It is clear that with the development of such advanced technologies, 
substantial improvements can be made in system size and weight. 

Concepts currently under consideration are for applications in 
low earth orbit (LEO) . These include solar powered Brayton and 
Stirling cycles operation at peak cycle temperatures of between 
1000 and 1400 K. To ensure continuous power during the eclipse 
portion of the LEO, thermal energy storage is accomplished through 
the use of a suitable phase change material (PCM) . With proper 
selection of materials the use of latent energy storage will provide 
relatively compact, low weight storage at a uniform high temperature. 
Both fluoride salts and metallic alloys are being evaluated for 
these applications. Typically the fluoride salts provide large 


heats of fusion with suitable containment chemical compatibility; 
they are limited by poor thermal transport properties and large 
volume changes. Metallic alloys provide outstanding thermal 
transport properties and restricted volume changes but do not have 
demonstrated containment materials. 

The work described herein addresses a portion of that development 
effort. The goal is first to examine a number of alternative 
receiver designs to determine their basic thermal performance. 
Designs providing large thermal gradients or large temperature 
swings during orbit are susceptible to early mechanical failure. 
In addition significant temperature variations in the working fluid 
have an adverse impact on the thermal efficiency of the cycle. By 
reviewing the thermal performance of each basic design, the relative 
merits of the basic concepts are to be compared. In addition, the 
effect of thermal enhancement and metal utilization provides a 
partial characterization of the performance improvements to be 
achieved by developing these technologies. The results are then 
thought to be applicable in determining the justification for 
development of new technologies. Finally the analysis includes a 
basic geometrical parametric study of each of the designs. Near 
the optimal design point it appears that there is generally a trade 
off between thermal performance and receiver weight. The results 
of this study will show how design parameters thermally affect 
overall system weight. 

1.1 Advanced Solar Dynamic Systems 

Initially the ASD program gave consideration to three basic 
power cycles: the Rankine, Brayton and Stirling engines. An 
organic Rankine cycle, using toluene as a working fluid, was 
considered but determined to be of limited potential for advanced 
applications. Organics tend to break down chemically at elevated 
temperatures, limiting the cycle to about 700 K. Since 
improvements in thermodynamic efficiency can most easily be 
achieved by increases in peak cycle temperature, consideration 
of the Rankine cycle then turned to either using liquid metals 
as working fluids or to finding an alternative stable, high 
temperature fluid. Preliminary design studies have shown that 
such systems are of significantly greater specific mass than 
their alternatives so that attention is currently concentrating 
on future development of the Brayton and 1'tirling cycles. 

1.1.1 Advanced Solar Dynamic Brayton Cycle 

A schematic of the Brayton cycle is shown in Figure 1. 
The cycle utilizes a mixture of xenon and helium as the working 
fluid. Earlier studies have indicated that a suitable 
combination of high and low molecular weight gases provides 
an optimum mixture with relatively high density and high thermal 
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conductivity [2]. The gas is compressed in a single stage 
mixed flow compressor, heated as it is passed successively 
through a recuperator and a receiver, expanded through a turbine 
and finally cooled by returning through the recuperator and 
passing through the heat exchanger. Energy is supplied by 
focusing solar energy from a parabolic collector into the 
receiver; waste heat is radiated into space. 



Figure 1. Advanced Solar Dynamic Brayton Cycle. 

1.1.2 Advanced Solar Dynamic Stirling Cycle 

The Stirling cycle is represented schematically in Figure 
2. The engine includes a piston and cylinder arrangement in 
which the working fluid moves alternately to chambers above 
and below the piston via a set of bypass ports in the cylinder 
walls. The upper chamber is thermally connected to the radiator. 
The piston is attached directly to a linear alternator to 
provide for a very simple system of electrical power generation. 

The Stirling engine represents a technological step beyond 
the Brayton engine. Whereas the Brayton cycle is widely used 
in terrestrial applications the Stirling engine is still in 
its infancy. The advantage lies primarily in its high specific 
power in larger generator sizes. Results of an earlier study 
[4] show that a 7 kWe Stirling system will weigh about 20% 
less than a comparable Brayton system, 40% less than a planer 
photovoltaic array, and 45% less than a photovoltaic array 
with a focused concentrator; at 35 kWe the advantage increases 




to about 25% below the Brayton cycle, 50% less than a photo- 
voltaic planer array, and 85% less than a photovoltaic array 
with a focused concentrator. NASA is currently testing a 
prototype Stirling engine to obtain additional data for 
evaluating its suitability for long term space applications. 



Figure 2. Advanced Solar Dynamic Stirling Cycle. 

1.2 Advanced Solar Dynamic Receiver Design 

Both the Brayton and Stirling cycles are adaptable to a wide 
range of high power applications. In low earth orbit provisions 
are made for uninterrupted power generation during the eclipse 
portion of the orbit by storing excess energy in periods of solar 
incidence as a sensible and/or latent heat. A receiver containing 
the storage media and designed for use with the Space Station is 
shown in Figure 3. On the right side of the receiver is an 
aperture placed at the focal point of the concentrator. Energy 
passing through the aperture is absorbed in the elements containing 
the PCM located on the inside perimeter of the receiver; a portion 
is stored for later use while the remainder is transferred to 
the working fluid flowing through a concentric passage inside 
the PCM. 


OReGINAL RAGfc IS 
OF POOR QUALITY 





Figure 3. Space Station Receiver Design. [3] 

Design studies [4] initiated to determine the overall weight 
distribution within a typical solar dynamic system have clearly 
established design targets for substantial performance 
improvements. Results, shown in Figure 4, indicate that the 
receiver would likely be the single heaviest component with about 
45% of the overall system weight. This is followed by the 
concentrator, with about 20% of the system weight, and the 
structure and radiator, each constituting about 10%. In an 
attempt to reduce overall system weight, it is natural that 
attention focus first on the receiver. In so doing it is also 
important to note that in most receiver designs, the PCM con- 
stitutes only a small fraction of the actual receiver weight. 
Receiver weight distributions for the base case design are shown 
in Table I. Here the PCM constitutes only 70 kg out of an overall 
receiver weight of 461 kg. This is a meer 15% and clearly suggests 
substantial opportunity for system weight savings through improved 
receiver design. 
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Figure 4. Weight Distribution for Space Station Solar Dynamic 

System. 


1.3 Project Scope. 

This study is designed to compare the basic thermal performance 
of five receiver concepts. The goal is to determine the efficiency 
of each in utilization of the PCM, the magnitude of the thermal 
gradients resulting, the extent of the resultant peak cycle 
temperature swings and their impact on overall cycle efficiency. 
Concepts were selected for analysis based partially on their 
respective merit, but also with the goal of including the widest 
possible range of distinct design approaches. Those selected 
include (a) a base case taken from the Space Station design 
effort, (b) a heat pipe configuration, (c) a pebble bed concept, 
(d) a Sanders Brayton cycle design and (e) a cavity heat pipe 
arrangement. The base case design is shown in Figure 3. The 
remaining concepts are shown in Figures 5 through 8 respectively. 











8 





Figure 6. Pebble Bed Receiver. [3] 

Each concept has been proposed by one or more outside con- 
tractors to NASA. However, the approach undertaken here is not 
to attempt to compare specific designs. These designs have been 
used to identify basic configurations. Each configuration has 
then been normalized to utilize the same group of PCMs and 
evaluated under a similar set of operation conditions irrespective 
of those originally proposed. The results are not seen as a 
basis for comparison of the respective proposals but instead are 
intended as a thermal comparison of the overall configurations. 

In addition to considerations of basic receiver configuration, 
optimal selection of those dimensional parameters which affect 
thermal performance will be of concern. Specifically this study 
will include a parametric study of each of the basic receiver 
designs to establish the range of performance achievable with 
each. While no direct trade-off will be performed between 
improvements in thermal performance and possible adverse effects 
on receiver weight, it is envisioned that the results might be 
useful in such a comparison. 
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Figure 7. Cavity Heat Pipe Receiver. [5] 

The third consideration in this study i ' that of PCM selection. 
Thermodynamic considerations would normally lead to a design in 
which all heating would occur at a temperature as high as practical 
below that imposed by the physical characteristics of the high 
temperature structure. This has led to the search for a material 
having a melting point near that imposed by the strength lim- 
itations of existing high temperature alloys and simultaneously 
possessing a high heat of fusion. 

A preliminary list of candidate phase change materials 
indicates numerous possibilities. The choice is complicated by 
the fact that many of these materials do not have well defined 
thermo-physical properties [6]. The characteristics sought here 
are not generally duplicated in other applications so that there 
has not heretofore been an intensive effort to establish prop- 
erties. Moreover the high temperature range of interest makes 
measurement particularly difficult. At such high temperatures 
materials are frequently quite corrosive making containment itself 
a problem. This analysis has been undertaken to determine the 
relative merits in utilization of either thermal conductivity 






enhanced salts or of metallic PCMs. By including direct eval- 
uations of the thermal performance of the respective receiver 
concepts incorporating such alternative materials is possible to 
ascertain the potential for improved performance from each. 
Ultimately such results may offer insight or lend weight to the 
decision to pursue development of the more promising technologies. 



Figure 8. Brayton Cycle Cavity Receiver. [5] 
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2 Analysis Approach 
2.1 Basic Cycles 

Performance parameters for the cycle have been established 
from typical Brayton Cycle conditions. The selection of the 
Brayton cycle appears natural in that the receivers selected for 
comparison were, with the exception of the Cavity Heat Pipe and 
the Stirling Heat Pipe designs, originally designed for that 
cycle. It is not thought that this assumption will significantly 
detract from the generality of the results. There appears little 
reason to believe that operating conditions would change dra- 
matically with either cycle so that results are thought to be 
generally applicable. 

The basic cycle performance characteristics are shown in 
Table II. The base design for comparison is a Brayton cycle 
engine producing 7 kw e from a continuous 25.9 kw^ input. This 
corresponds to an overall cycle efficiency of 27%. Variations 
in the turbine inlet temperature will affect both cycle efficiency 
and power output; these effects may probably be best handled by 
controlling the working fluid flow rate or bypassing a portion 
of the flow around the receiver. In any case these control 
mechanisms are not used for purposes of this analysis and both 
receiver inlet temperature and flow rate are maintained at a 
constant level. The coolant composition, inlet temperature and 
flow rate are taken as the nominal values used in the base design. 
This selection might be considered to favor the base case but in 
preliminary studies, designs have not shown a significant sen- 
sitivity to the selection. In any case it was desired to compare 
the respective designs under similar operating conditions so that 
utilization of base case parameters appeared more reasonable than 
selecting conditions proposed with one of the advanced designs. 

TABLE II. BASIC CYCLE PERFORMANCE CHARACTERISTICS 


Parameter 

Specification 

Engine Power 

7 kw e 

Engine Thermal Input 

2 5.9 kw-(- 

Receiver Losses 

3.1 kw t 

Solar Incidence Period 

54 min. 

Eclipse Period 

36 min. 

Working Fluid 

71% He 29% Xe 

Coolant Flow Rate 

0.25 kg/s 

Coolant Inlet Tempera- 

800 K 

ture 





2 . 2 Concentrators 

The performance of any receiver in which the PCM is directly 
irradiated by the incident solar flux will be strongly a function 
of the incident flux distribution. In the study reported here 
all flux distribution have been determined from the CAV2 [7] 
computer program. The input parameters used in each of the cases 
are as given in Table III. The receiver diameter and depth were 
selected to accommodate the PCM elements in each case without 
additional clearance. The elements were located at a distance 
from the aperture so that the incident flux would begin to become 
significant at the first PCM cannister and was extended to exactly 
the element length. The inside diameter of the receiver was set 
at the outer diameter formed by the PCM elements. 

TABLE III. BASIC CONCENTRATOR CHARACTERISTICS 


Parameter 

Specification 

Dish diameter 

3.59 m 

Rim angle 

45° 

Slope error 

1% 

Aperture radius 

0.091 m 

Aperture offset 

0 m 

Dish misallignment 

0 mr 



Figure 9. Base Case Incident Solar Flux Distribution. 
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The resultant flux distribution calculated for the base case 
configuration is shown in Figure 9. The distribution is shown 
for receiver arrangements with L/D ratios varing from 1 to 2 as 
used in the parametric evaluation. Clearly the flux is quite 
high near the aperture end of the receiver leading to difficulties 
in obtaining uniform PCM utilization. The problem is somewhat 
more severe with the higher L/D ratios indicating a higher incident 
flux to compensate for the lower reflectance from the bottom of 
the receiver. Receivers utilizing a truncated cone configuration 
can provide a more uniform incident flux distribution and a more 
appropriate PCM allocation. 

The shape of the incident flux curves, while relatively 
simple in form, are somewhat difficult to represent numer- 
ically without introducing certain round off error. In order 
to smooth the function it has been decided to integrate the 
area under the curve allowing representation of the cumulative 
heat input as a function of axial position. The integration 
has been performed using standard numerical techniques and 
resultant curve is shown in Figure 10. Also shown is a third 
order curve fit of the integral used in the receiver program. 
The technique allows for relatively accurate representation of 
the CAV2 results but also provides a somewhat clearer picture 
of the tendency for poor flux distribution. It may be noted 
that 48% of the incident energy enters the first 20% of the 
receiver length. 



Figure 10. Receiver Energy Distribution. 
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2 . 3 PCM Enhancement 

While fluoride salts have been considered as prime candidates 
for selection as the phase change material, they are characterized 
by a low thermal conductivity. In order to achieve a high level 
of thermal utilization and low temperature gradients designers 
have relied on a variety of thermal enhancement schemes. Proposed 
schemes have included high conductivity metal fins extending 
through the PCM, wire meshes within the PCM, metallic foams in 
which the PCM can be placed and a graphite matrix in which 
individual PCM pellets are to be arrayed. Manufacturing of 
certain of these concepts represents an unproven technology. 
Also unproven is the long term integrity of the concepts during 
normal orbital cycles. Nevertheless there remains the need to 
develop techniques of improving the thermal performance of such 
materials and it would seem prudent at this time to evaluate the 
impact of specific material enhancement schemes so as to perhaps 
limit development to only the most promising concepts. 

The approach used here has been to develop a model based 
largely around a "random" metal mesh. It is probably easiest to 
envision such a mesh as being constructed from layers of wires 
lying in a plane. Stacked on the first layer will be a second 
layer of wires arranged so that the wires in adjacent lamina 
appear to form a square grid. The model is designed around a 
large number of such planes. If one of these laminas is viewed 
from the edge so as to be seen from the end of the wires, the 
wires and PCM may be considered to act to transfer heat in 
parallel. The thermal resistance associated with this plane is: 

R 2 

11 ( 1 ~ £ ) ' R PCM + ■ k F. :V // 


Adjacent planes, when viewed from the same perspective, will 
appear as layers of PCM acting in series with the wire segments. 
The effective thermal resistance of these laminae are given as: 




P-2_r 
0 .5 • k pr i 


P 


The composite system will appear as a large number of such planes 
all acting in parallel with one another: 


/?. R,*R S 

When this system is viewed from a direction normal to both layers 
of wires, each lamina will appear as a cylinder in parallel with 
PCM. The equivalent resistance of this arrangement is given as: 
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J k pc At * f\ H ( P ~ 2 • r ) ■ k PCM 

l 0.5 



c ?C.U f,v« (P-2r)A -PCM 


It is assumed that these two orientations are randomly dispersed 
within the elements so that the overall effect is to produce an 
isotropic effective conductivity. It is further assumed that 
these resistances will combine in the same manner as structural 
properties in a composite material. From random orientations of 
fibers in a composite materials: 


EFFECTIVE 


8 R , 8 R 


This approach has be utilized to evaluate the effect of adding 
the various enhancement materials. The results are shown in 
Table IV. In the listing, the density is given as the density 
of PCM per unit volume. Since voids are assumed to remain 
distributed within the solid phase the same value is reported 
for liquid and solid. Likewise, the specific heat reported is 
based solely on the PCM so that the enhancement material does 
not add to sensible energy storage capacity. The thermal 
conductivity values include the effects of distributed voids 
within the solid and may be taken as the effective enhanced 
conductivity . 

Table IV. Thermal Transport Properties of Enhanced PCM 


Material 

Density 

kJ/m 3 

Specific 
Heat 
kj/kg K 
(Sol/liq) 

Thermal 
Conductivity 
kw/m K 
(Sol/liq) 

Heat of 
Fusion 
kJ/kg 

Energy 

Density 

gj/m 3 

LiF 22 CaF 2 

2097 . 0 

1.84/1.97 

. 00413/. 00173 

753 

1.597 

10% BN 

1887.3 

1.84/1.97 

.00492/. 00254 

753 

1.421 

18% BN 

1719.4 

1.84/1.97 

. 00554/ . 00320 

753 

1.295 

10% Grap. 

1887 . 3 

1.84/1.97 

.00643/. 00387 

753 

1.421 

18% Grap. 

1719.4 

1.84/1.97 

. 00825/. 00559 

753 

1.295 

Ge 

5010.2 

0.38/0.39 

.0405/. 0424 

481 

2.410 


2.4 Receiver Concepts 

2.4.1 Base Receiver Concept 

The base case receiver design is essentially similar to 
that shown in Figure 3 for the space station. In this design 













the PCM is contained in individual canisters. These canisters 
are constructed as a concentric annulus with each end closed. 
Several of these individual canisters are then stacked by 
placing them over a central coolant tube. The coolant tube 
with the canisters forms a single PCM element. Several elements 
are then placed inside the inner circumference of the receiver 
shell with equal azimuthal spacing. The coolant tubes are 
connected by separate manifolds at each end of the receiver 
to form the inlet and outlet gas passages. Dimensions on the 
base case receiver are shown in Table IV. 


TABLE V. BASE CASE RECEIVER CHARACTERISTICS 


Parameter 

Specification 

Tube Length, m 

1.8 

Number of Tubes 

75 

Sidewall/PCM Thickness 

12.5% 

Coolant Tube ID 

0.00841 m 

Canister ID 

0.01095 m 

Canister OD 

0.015 m 

Canister Wall Thickness 

0.00071 m 

Element P/D 

1.66 

PCM 

LiF 22 CaF 2 

PCM Mass 

96 kg 


2.4.2 Heat Pipe Concept 

The heat pipe receiver design is shown in Figure 5. The 
receiver is constructed in three axial sections. The first 
of these sections is located nearest the receiver aperture and 
includes no PCM canisters. Solar energy entering through the 
aperture strikes the central tubes directly. Inside the tube 
is a heat pipe working fluid. The evaporated fluid then 
transports the energy into the remaining two sections. In the 
central section the tube is surrounded by PCM canisters, with 
a configuration much like those in the base case design. A 
portion of the energy from the heat pipe fluid is conducted 
into the PCM and stored as latent energy during periods of 
solar incidence. The remaining energy from the heat pipe is 
transported to the last section of the receiver where a heat 
exchanger couples the heat pipe to the engine. During the 
solar eclipse the heat pipe drops in pressure and temperature 
and begins to remove energy from the PCM and transports it to 
the heat exchanger in the last section of the receiver. 
Dimensions on the base case receiver are shown in Table V. 




TABLE VI. HEAT PIPE RECEIVER CHARACTERISTICS 


Parameter 

Specification 

Tube Length, m 
Number of Tubes 
Sidewall/PCM Thickness 
Coolant Tube ID 
Canister ID 
Canister OD 

Canister Wall Thickness 
PCM 

PCM Mass 

0.7116 m 
20 
12.5% 
0.0482 m 
0.0522 m 
0.0820 m 
0.00127 m 
LiF 22 CaF 2 
91.5 kg 


2.4.3 Pebble Bed Concept 

The pebble bed receiver is shown in Figure 6. In this 
case the PCM is stored in individual capsules, each coated 
with a suitable high temperature material. The capsules are 
contained within the coolant passage. The outside surface of 
the coolant tubes is directly exposed to the solar flux. Energy 
is then conducted across the tube wall where a portion is 
convected directly to the cycle working fluid and the remainder 
is conducted to those solid particles adjacent to the wall. 
Energy is transported to the interior elements by a combination 
of conduction and convection through the working fluid. 
Dimensions for the pebble bed concept are shown in Table VI. 

TABLE VI. PEBBLE BED RECEIVER CHARACTERISTICS 


Parameter 

Specification 

Tube Length, m 
Number of Tubes 
Packing Fraction 
Coolant Tube ID 
Canister Wall Thickness 
PCM 

PCM Mass 
PCM Sphere OD 

0.91 m 
15 
70.0% 
0.82 m 
0.00127 m 
LiF 22 CaF 2 
96 kg 
0.014 m 


2.4.4 Cavity Heat Pipe 

The system under consideration will consist of a cavity 
heat pipe, shown in the Figure 7. This design is intended 
for use in conjunction with the Stirling engine, but it is 
believed that a similar design could be developed for the 




Brayton engine with the inclusion of a suitable heat exchanger 
within the receiver. In this design the receiver is configured 
as a right cylinder. The solar incidence from a parabolic 
collector arrives at the receiver and passes through an aperture 
at the bottom. Inside the aperture the flux strikes the 
evaporator section of the cavity. The vapor then transfers 
heat to condensing surfaces on both the thermal storage material 
(TSM) and the engine high temperature heat exchanger. 

The thermal storage material is contained in a number of 
canisters arranged around the periphery of the receiver. The 
side walls of these structures are tapered inward toward the 
apex, so that as the canisters are placed side by side they 
bend into a closed cylinder configuration. This cylinder is 
then banded together and placed inside the receiver. The apex 
of the canisters are covered with a condenser/ evaporator 
material, generally a wire mesh. During the period of solar 
incidence, this area receives a continuous heat input, gradually 
melting the TSM. During the eclipse portion of the orbit, 
energy is removed from the TSM by evaporating the heat pipe 
fluid from this area. 

The shape of individual canisters has been established as 
a method of thermal enhancement. In effect, if one were to 
consider the region between roofs, this region acts as a 
triangular shaped fin with the heat pipe fluid replacing the 
fin material. The arrangement provides for fins of extremely 
high effective conductivity and of relatively low equivalent 
mass. Increasing the number of fins will necessarily imply 
an increased number of canisters, so that careful trade-offs 
are necessary between thermal enhancement and system weight. 

The high apex shape may also serve to alleviate the problem 
of void management, but this has not yet been confirmed. 
Certainly the narrowing walls are reminiscent of geometries 
proposed for fuel management within zero ?g environments. 
However Marangoni forces may well counter these effects during 
the heating phase and a careful analysis will be required to 
ascertain the final deposition. 

From a thermal aspect, the problem is not so severe as in 
those designs which allow direct solar incidence onto the 
canisters. In those cases the heat flux at any point in the 
receiver was fixed by the collector and receiver geometries. 
If there were inadequate material to absorb the incident 
radiation, surface temperatures would rise and a thermally 
induced failure might occur. In the cavity heat pipe design 
it is not the flux on the canister surface which is fixed, but 
rather the temperature. The temperature throughout the cavity 
will be nearly uniform and adequate to drive the energy into 



the TSM, neither more nor less. At any point where low thermal 
storage capacity might occur, the energy input would simply 
decrease locally and the corresponding energy would be stored 
elsewhere. Only in the case where the void covered a significant 
fraction of the overall heat transfer surface would a problem 
arise. In this analysis it is assumed that massive void 
coverage of the heated surface will be avoided. 

TABLE VII. CAVITY HEAT PIPE CHARACTERISTICS 


Parameter 

Specification 

Tube Length, m 
Number of Tubes 
Packing Fraction 
Coolant Tube ID 
Canister Wall Thickness 
PCM 

PCM Mass 

0.91m 
15 
70.0% 
0.82 m 
0.00127 m 
LiF 22 CaF 2 
96 kg 


2.4.5 Brayton Cycle Cavity Receiver 

The Brayton cycle cavity receiver concept is shown in 
Figure 8. A major feature of this particular design are a 
number of innovative concepts to decrease reciever weight. 
Most notable is the method used for PCM containment. The 
overall configuration of the receiver is that of a truncated 
cone with corregations running from the base toward the apex. 
The PCM is held into the troughs of the corregations by capillary 
action while being contained within the receiver by a suitable 
window across the receiver aperture. A second smooth conical 
section is situated outside the first, with the interstitial 
space then serving as a passage for coolant. This particular 
arrangement is interesting in that it provides for a relatively 
light structure and includes a relatively limited number of 
fabrication steps. 

Solar radiation which enters the receiver inpinges directly 
on the PCM. Since the PCM is normally highly transparent to 
the incident radiation, suitable addatives are to be included 
to ensure adequate absorption. While this approach is as yet 
unproven, it is interesting from a heat transfer viewpoint in 
that it allows the designer a degree of control over the 
absorbed energy distribution. An indirect consequence of this 
is that high element P/D ratios are not required to achieve 
effective distribution of the stored energy. The net result 
is an extremely compact receiver with assocated weight savings 
in structure and insulation. 



This concept is the only one of the designs considered 
which utilizes the truncated cone geometry to distribute PCM 
with higher fractions near the aperture and a decreasing 
fraction at the more distant locations. By situating the PCM 
where it can be most effectively utilized a somewhat decreased 
mass is required. 

TABLE IX. BRAYTON CYCLE RECEIVER CHARACTERISTICS 


Parameter 

Specification 

Corrugation Length 
Number of Corrugations 
Height of Corrugation 
Corrugation Wall Thickness 
Base Radius of Truncated Cone 
Top Radius of Truncated Cone 
PCM 

PCM Mass 

0.654 m 
100 

.064127 m 
. 001 m 
.386 m 
.154 m 
Lif 22 Caf 
96 Kg 


2.5 Hybrid Elements 

Within the parametric studies performed on each of the receiver 
concepts, four PCM concepts have been evaluated. These have 
included (a) a fluoride salt, (b) a fluoride salt enhanced with 
Boron Nitride, (c) a fluoride salt enhanced with ATJ graphite 
and (d) a metallic PCM. This grouping, while considering only 
a limited number of alternatives, is thought to bracket the range 
of thermal transport properties of most materials under con- 
sideration for Advanced Solar Dynamic applications at this time. 

As a part of this study several additional materials were 
considered and, while not included in the parametric evaluation, 
were evaluated for thermal properties. In this case NaF was 
selected as the fluoride salt with Germanium being selected as 
the metallic alternative. Enhanced combinations each included 
NaF with the following enhancement schemes: (a) 12% of volume 
consisting of Niobium fins, (b) 12% Nickel fins, (c) 18% niobium 
wire mesh, (d) 18% nickel wire mesh and (e) 18% pyrolytic graphite 
mesh. 

An alternative arrangement considered was to encapsulate 
Germanium into spheres and to fill the interstices in the packing 
with either a conductivity enhancement material, i.e. graphite, 
or with a nonencapsulated PCM, i.e. NaF. Using uniform sized 
spheres it is possible to achieve packings of spheres of about 
60%. By using 2 significantly different sphere sizes packing 
densities of up to about 80% are achievable. 




A final arrangement considered included packings of uniform 
sized spheres of encapsulated Germanium with pyrolytic graphite 
fibers traveling along the major axis of the packing to enhance 
thermal conductivity. The interstitial spaces were then filled 
with NaF to completely utilized the available volume. 

The results of an analysis of each of these schemes to 
determine the effective thermal transport properties is shown in 
Table X. In each of these comparisons the most important thermal 
properties are the effective thermal conductivity and the energy 
density. The effective conductivity is of course that hypothetical 
thermal conductivity which would produce the same temperature 
difference in a large sample as would the nonhomogeneous material 
under consideration. It may be viewed as indicative of the 
temperature gradients which each sample would produce in an 
operating system. The energy density is the latent energy stored 
per unit volume by the nonhomogeneous element. It may be viewed 
as an indication of the required size of receiver which might 
incorporate that particular PCM. 

What is not indicated within this comparison is the ease of 
development of each of the respective concepts. Whereas Germanium 
appears to be the choice candidate on the basis of energy density 
and thermal conductivity, problems of containment may suggest an 
alternative selection. However almost any other selection of 
Germanium in combination with other materials considered here 
tends to result in an improved energy density. The potential 
combination with NaF, if it could be fabricated, certainly 
indicates a potential for compact design. However this particular 
combination results in an unusually low thermal conductivity and 
it not thought to merit further consideration. 

The enhancement of NaF with either of the two graphites, ATJ 
or pyrolytic graphite, will serve to significantly improve the 
effective thermal conductivity. The use of pyrolytic graphite 
offers a significant advantage over nickel additions. However 
any addition of a nonmelting material will necessarily decrease 
the energy density. The choice of a fluoride salt with any of 
the thermal enhancement additions is a low risk choice in terms 
of technological development while offering satisfactory thermal 
conductivity improvement. It will, of course, require a relatively 
large and heavy receiver. 



Table X. Thermal Transport Properties of Advanced Enhancement 
PCM 


Material 

Density 

kJ/m 3 

Specific 
Heat 
kJ/kg K 
(Sol/liq) 

Thermal 
Conductivity 
kW/m K 
(Sol/liq) 

Heat of 
Fusion 
kJ/kg 

Energy 

Density 

gJ/m 3 

NaF 

1948.6 

1.12/1.24 

.003/. 00013 

798 

1.537 

G© 

5010.2 

. 381/. 396 

.0405/. 0124 

481 

2 .401 

Fins : 






NaF/Nio. 

2744.4 

1.01/ 

.0087/. 0069 

493 

1.360 

NaF/Ni. 

2744.4 

1.01 

.0105/. 0081 

493 

1.360 

Mesh : 






NaF/Nio. 

3142.3 

.964/1.06 

.0044/. 0035 

401 

1.261 

NaF/Ni. 

3142 . 3 

.967/1.07 

.0052/. 0042 

401 

1.261 

NaF/ Pyro. 




401 

1.261 

Matrix: 






. 4NaF/ . 6Ge 

3785.6 

.680/. 730 

.0045/. 0021 

545 

2 . 061 

. 2NaF/ . 8Ge 

4397.9 

.530/. 560 

.0053/. 0026 

508 

2 .235 

. 6Ge/Pyro. 

3685.3 

.90/. 91 

.0559/. 0577 

392 

1.446 

. 8Ge/Pyro . 

4347 . 7 

. 640/ . 652 

. 0340/. 0365 

443 

1.926 

. 8Ge/ATJ 

4347 . 7 

.640/. 652 

.0188/. 0203 

443 

1.926 

Hybrid: 






NaF/Ge/Pyro 

3864.9 

.499/. 516 

. 0204/. 0179 

455 

1.757 


By a process of elimination it would appear that the 
aggressive, high payback choice would be a combination of Germanium 
encapsulated or contained in graphite. If high density, i.e. 
80%, packing can be achieved this approach combines both unusually 
high thermal conductivity and high energy density. 

2.6 Model 

In this study effort, a number of semi-two dimensional, 
transient, finite difference code were developed to evaluate each 
of the receiver systems. The codes were then used to perform 
several pertinent parametric studies. These studies were directed 
not only as establishing the performance of the proposed systems, 
but also toward optimizing such performance and establishing 
directions of future development. 









While describing the model, reference will be made to the 
base line configuration. In certain of the other concepts this 
approach has been modified to represent the specific arrangements, 
but in general the overall methods remain uniform throughout. 
The approach has been to utilize a two dimensional, transient 
analysis with an implicit formulation. All conduction is assumed 
to occur in the axial (R-Theta or X-Y, depending on the coordinate 
system selected) plane. A series of such axial elements are 
stacked between the coolant inlet and outlet. For each axial 
element the solar input is determined based on the results of 
the CAV2 output. Heat input into the coolant is determined based 
on the calculated convective coefficient, coolant inlet tem- 
perature and the average wall temperature. Successive axial 
elements are connected thermally through the change in coolant 
temperature as it flows downstream. This approach utilizing the 
implicit calculational scheme provides for stable solutions for 
relatively long time steps and therefore is thought to be most 
applicable to computations involving medium to small computers. 

The two dimensional finite difference program is designed to 
describe the system shown in Figure" 11. Symmetry is assumed so 
that only one half of the element is modeled. Within the 
semi-annulus an arbitrary number of elements is used in the R 
and Theta directions. Generally this will require annular 
segimental elements where the dR segiments are uniformly spaced 
and, in general, unequal to dTheta. Within the metal containment 
representing the first and last radial elements, a somewhat 
different radial length is used corresponding to the actual 
thickness of the containment wall. 

All radiant heat transfer is assumed to occur through the 
outside circumference. In general this flux will be nonuniformly 
distributed azimuthally. In the model the total incident flux 
is constituted from a uniformly distributed portion representing 
the indirect radiation and a frontal portion representing a direct 
radiation component. Because of symmetry both the 0 and 180 
degree walls are insulated. All convective losses are assumed 
to occur through the inside circumference. Convective coeffi- 
cients are determined from conventional pipe flow correlations 
and assume a developing thermal boundary layer based on the 
distance from the inlet to the first element. 

There are, within the literature, a number of methods of 
presenting the energy equation for solution of problems involving 
phase change. The method employed here is not taken directly 
from any of these but is a logical extension of the methods 
proposed by Pantankar [8] for the general solution of conduction 
problems. Emphasis is placed on writing relatively short, simple 
programs using methods that are broadly understood within the 
heat transfer community. While it might be argued that certain 
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of the advanced procedures provide for improved accuracy or 
stability, the advantages of presenting a scheme generally 
understood within the potential range of users or reviewers and 
ease of technical support are considered to be compelling. 



Figure 11. Finite Difference Element Arrangement. 

2.6.1 Physical model. 

The physical model is somewhat complicated by a number of 
factors including zero g effects, inadequate physical property 
data and the need for generality where various thermal 
enhancement devices might be considered. Particular caution 
must be employed in such choices in that the selection of 
differing models for alternative conditions may, in fact, 
overshadow all other parameters under consideration. This is 
particularly true in terms of the void model and while alternate 
enhancement schemes may actually justify differing models a 
uniform model has been incorporated here in an attempt at 
generality. 

2.6.2 Void model. 

Considerable uncertainty exists as to void location during 
the heating/cooling cycle. A characteristic of the fluoride 
salts currently being considered is that they tend to undergo 
significant density changes with phase change. Normally this 
will involve a density increase of between 10 and 25% as the 
material solidifies. Thus upon solidification voids will tend 




to form so that the question remains as to how such voids 
react. If they fail to coalesce then such voids will likely 
remain dispersed within the solid phase. If they do coalesce 
then the question remains as to how such voids migrate and to 
what position. 

At this time an experiment is planned for a space flight 
test in an attempt to gain actual data on void migration. 
Prior to acquiring such data any model will be, at best, only 
speculative. In those cases where a high conductivity foam 
is to be used for enhancement of the thermal conductivity, it 
might reasonably be anticipated that migration will be 
restricted and that voids will remain within the cells where 
they are formed. It is thought possible that viscous effects 
might produce similar distributions in the more general case. 

In the absence of experimental data it is assumed that all 
voids form microspheres which are uniformly distributed within 
the solid. The remaining material will then exhibit no 
convective component to heat transfer and will retain a uniform 
effective density throughout. Note that the thermal conduc- 
tivity will be affected. In this case the granular material 
model proposed by Meredith and Tobias [9] is utilized to predict 
the effective thermal conductivity of the overall system. For 
the case in which radiation is neglected, the voids may be 
assumed to have zero thermal conductivity and the equation 
reduces to the form: 

k. ff 2- ( 1 - c)-0.6 135- c 2333 - 1 . 6 • (•■ 3333 

A: tsai ~ 2 + t + 0.6 135- e 2 333 - 0.69/5 -e 3 333 


where c = the void fraction. 

An alternate approach is used in the case of the cavity 
receiver. In this case if voids are assumed to coalesce they 
may be assumed to migrate to the liquid-cavity interface and 
to leave a 100% dense solid phase. This is in contrast to 
those enclosed PCM systems where the final location of void 
migration remains unanswered. To remain consistent with the 
analysis of the original contractor proposing this concept it 
was decided to assume such migration in this one case. This 
approach will tend to result in higher calculated thermal 
conductivities between the melt front and the gas coolant. 
The effect will clearly be most significant when the PCM is 
largely solidified. This occurs early in the heating cycle 
and late in the cooling cycle. It should be noted that the 
time early in the heating cycle does not correspond to a time 
of temperature extremes so that the consequences on engine 
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performance and thermal stresses are not thought to be sig- 
nificant. However, the period late in the cooling cycle is a 
time of temperature extremes. The minimum predicted turbine 
inlet temperatures may be somewhat nonconservative. Thermal 
stresses at this same point in time are not thought to be as 
severe as those occurring late in the heating cycle so that 
these results remain uneffected. 

Computationally the void migration is accomplished by 
allowing a mass interchange between any element undergoing a 
phase change and the outermost element within the same radial 
segment. At each time step the movement of the melt front is 
calculated within the element undergoing phase change. The 
change in mass is then calculated. It is further assumed that 
the transferred mass leaves and enters the respective elements 
at the melting temperature as a method of ensuring that 
appropriate energy balances are maintained. 

2.6.3 Energy equation. 

The general transient energy equation for a stationary 
system may be written in terms of sensible and latent energy 
accumulation and conductive transport as: 

cLT 7 ,_ /7 . 
pc„ — — — + ph fl w k ■ n 


where ?/ = melt front velocity. 


In finite difference form this equation may be written as: 


Ax Ay Ax 


1 = 2 k±±(T„ 
Ay 


T,„y" *2k±±{T R + T,y" 

Ay 


+ — T T + p h ,, ( 
AT ' 


The major feature in descretizing this formula is in the 
treatment of the latent accumulation term. The resultant 
equation is nonlinear so that a trial and error procedure is 
required. Since it is anticipated that the temperature of any 
element with a partial phase change should be near the melting 
temperature, the velocity may be adjusted in each iteration 
so as to cause the associated node temperature to vanish. 


v = 


l-i c „ ( 7 7 Ine [( ) 


where (.i = a convergence factor. 



In the cases reviewed here the convergence factor could be 
taken as between 0.25 and 2.0. The velocity term obtained in 
this manner is cumulative. That is, it is underevaluated in 
each iteration so that summing the value for subsequent 
iterations will eventually lead to a convergent value where 
the local temperature approaches the melting value. A second 
convergence criteria is also included based on an energy balance 
for the overall axial segiment. In general the goal is to 
seek energy convergence to within . 1% . 

2.6.4 Computer Program. 

Selection of an implicit scheme to solve the nodal equations 
is designed to enhance numerical stability and to permit the 
use of reasonably long time increments. This is viewed as a 
necessity for such programs when utilizing relatively small 
computational facilities. A major advantage to such an approach 
is that the resultant program is highly portable, being designed 
primarily for PC usage. 

In order to simplify the programming and to reduce the 
number of parameters being stored a hybrid computational scheme 
is employed which included features of both Gaussian elimination 
and Gauss Seidel. Considering the figure of the finite 
difference scheme, the approach is to successively solve for 
the temperatures along each of the vertical columns, using a 
one dimensional (tridiagonal) method. While evaluating the 
temperatures along any one column, temperatures along any other 
column are assumed to be known based on the last iteration. 
Utilizing successive sweeps, this method provides for a rapid 
convergence while limiting memory requirements. 
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3 Parametric Study 

In the parametric studies which follow, the effects of the basic 
geometric parameters influencing overall thermal performance of the 
canisters are examined. These include the number of canisters, 
element geometry or spacing and the use of alternative thermal 
enhancement materials. During the transient, temperatures are 
evaluated at two times, 54 and 90 minutes. These times correspond 
to the end of the solar incidence, where maximum temperatures are 
achieved, and to the end of the solar eclipse, where minimum 
temperatures are encountered. The implications of such temperatures, 
while largely qualitative, directly influence both receiver life 
and cycle efficiency. 

One important aspect of the study is to underline the differences 
in thermal performance between metals and salts as PCM's. Because 
of the inherently low thermal conductivity of the salts, their 
parametric range has been extended with a view of establishing 
potential for improved performance with improved geometry. In the 
case of metals this same extension has appeared unnecessary. 

3.1 Base Line Receiver. 

The base line receiver utilizes convective transport directly 
to the working fluid rather than utilizing a heat pipe as the 
other designs considered here. The nature of the design constrains 
the parametric study to some extent so that certain combinations 
of enhancement and configuration perturbation are not possible. 
This is not a limitation of the design, but rather of the parameters 
selected. In order to take advantage of the full range of 
enhancement the system would require a different inlet temperature 
and/or flow rate. The values chosen for the study basically 
correspond to those initially proposed in the design. This 
decision has meant that considerations oii enhancement are thus 
addressed on a more limited scale. 

It should also be noted that any enhancement effect due to 
the canister side walls is neglected. If foam enhancement is 
added it is assumed that the foam will provide for the void 
management so that a large number of individual canisters is not 
required. Since exclusion of the canister walls is thought 
necessary for evaluation of enhancement effects it has also been 
omitted from the base case to fully indicate relative improvements. 

The turbine inlet temperature for the base line case during 
normal orbit is shown in Figure 12. Thermouynamically the outlet 
temperature must average 1000 K (LiF 22 CaF 2 melts at 1039 K) 
and variations around this value are due to alternate melting 
and sensible heating or solidifying and sensible cooling of the 
PCM. The temperature is seen to increase during the solar 



incidence reaching a value about 30 degrees above the average 
and cooling rapidly during eclipse to about 30 degrees below the 
average . 



Figure 12. Base Line Orbital Turbine Inlet Temperature. 

The effects of thermal enhancement are seen in Figure 13. 
Note that the plotted values do not extend to include the 18% 
enhancement value. A projection of the partial low temperature 
lines would indicate that this case would fall above the average 
gas temperature violating the first law. However it is seen 
that over the effective range the lower temperature is increased 
an average of 3.4° and 4.2° for each percent addition of BN and 
pyrolytic graphite, respectively and that the enhanced materials 
can be utilized effectively to control temperature swings. Full 
utilization of the enhancement would require either an increase 
in coolant flow rate or reoptimization of the inlet temperature. 

Increasing the number of elements serves to reduce the radius 
of the individual canisters and subsequently reduces the thermal 
resistance within each. As indicated in Figure 14 the effect is 
fairly significant up to between 60 and 75 elements producing 
about a 3° improvement per element. In view of the weight 
advantage of reducing the number of elements it would appear that 
the optimal condition would occur at a number near 60 or 65 
elements. 
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distribution was assumed to remain constant with each of the 
configurations. In fact as the L/D increases the heat flux will 
tend to become somewhat more maldistributed. In this case the 
improvements seen at high L/D ratios will be somewhat reduced. 



Figure 15. Effects of Receiver L/D on Turbine Inlet Tempera- 
ture. 

The effect of increasing the element circumferential pitch 
to diameter ratio would be to decrease the azimuthal variation 
in the incident heat flux. The improvement in the utilization 
of the PCM will then lead to an improvement in the uniformity of 
the inlet turbine gas temperature. This effect is seen in Figure 
16. It is noted that the reduction in variation appears nearly 
linear over the P/D range of 1.25 to 2.0 with a corresponding 
improvement of 30°. 

In Figures 17 through 19 the maximum and minimum canister 
temperatures are shown at 54 minutes in the first section of the 
element. These values are obtained from the same runs in which 
the turbine inlet temperatures were calculated. What will be 
generally noted is that the PCM temperatures show much the same 
trends but are much more sensitive to the geometrical parameters. 
Where poor PCM utilization is encountered there is a shift from 
storing energy in its latent form to storage as sensible heat. 
The inevitable results include large variations in turbine inlet 
temperature accompanied by high peak PCM temperatures. 

The PCM temperature differential as a result of increasing 
the element number is shown in Figure 17. The design is subjected 
to considerable peaking of the incident heat flux and shows 
considerable sensible heating until the number of elements nears 
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The effect of element pitch to diameter ratio on PCM tem- 
peratures is shown for the base line receiver in Figure 19. Peak 
temperatures decrease from about 270° to about 125° above melting 
as the P/D increases from 1.25 to 2. Canister temperature 
differences decrease from about 40° to about 20° over the same 
range. 

3.2 Heat Pipe Receiver. 



Figure 20. Effects of PCM Selection on Heat Pipe Vapor 

Temperature. 

In Figure 20 the heat pipe vapor temperature is indicated 
through the steady orbit. It can be seen that temperatures change 
markedly between eclipse and solar incidence. The heat pipe 
fluid must operate at above melting to store energy during the 
incidence period and below melting to remove energy during eclipse. 
As either period progresses the melt front advances through the 
PCM so that a continually increasing temperature difference is 
required to sustain the heat transfer rate. One can also see 
the effects of thermal enhancement of the PCM. Additions of 18% 
by volume of BN or ATJ graphite substantially reduce the indicated 
temperature variations with overall variations decreasing from 
about 110° to about 75° and 60° respectively. The use of a metal 
PCM is seen to virtually eliminate all such variations providing 
only about a 25° change. Note that virtually all of the variation 
in the last case occurs in the very last time step. This is an 
indication of excessive solidification and use of only sensible 
energy. With some adjustment of the operating conditions this 
drop can be eliminated reducing the overall variation to only 
about 5°. 
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The maximum and minimum heat pipe vapor temperatures cor- 
responding to the end of solar incidence (54 min.) and the end 
of solar eclipse (90 min.) as a function of the percent of thermal 
enhancement are shown in Figure 21. Also shown are the corre- 
sponding swings encountered with the metal PCM (Germanium) . 
Enhancement materials serve to improve the radial conduction but 
the effect is partially countered by displacing PCM. This produces 
an improved conductivity but requires that heat be conducted over 
a larger distance. While the effect is significant in either of 
the enhancement materials, both are substantially less effective 
than the use of metals. From the heat pipe vapor temperatures 
produced it would appear the between 8 and 15% of either enhancement 
material is required to produce a significant improvement. 



Figure 21. Effect of Enhancement Fraction on Heat Pipe Vapor 

Temperature. 

An increased number of elements has the effect, as seen in 
Figure 22, of decreasing the radial thickness of PCM in each 
canister. The result should be a decrease in heat pipe vapor 
temperature variation as the thermal resistance between the heat 
pipe fluid and the melt front is decreased. With LiF 22 CaF 2 
the variation in heat pipe vapor temperatures is seen to be 
substantially reduced as the number of elements is increased to 
between about 20 and 25. Beyond this point the impact of further 
additions appears limited. The effect can be achieved with far 
fewer elements with some degree of thermal enhancement. The 
metal PCM is shown only for the case of 15 elements. It would 
be expected that the effect of increasing the canister number 
for the Ge PCM would result in a similar curve for the salts. 





Figure 22. Effect of Element Number on Heat Pipe Vapor 

Temperatures . 



Figure 23. Effect of Receiver L/D Ratio on Heat Pipe Vapor 

Temperatures . 

As the receiver is made longer the radial dimension of the 
individual canisters is decreased. Thermally this results in a 
decrease in both the heat flux and the radial resistance through 
which energy must flow. As seen in Figure 23 the effect appears 
to be relatively large for ratios up to 1. For longer receivers 
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there appears to be only a small improvement in the variation. 
Again the effect of the solar flux profile would be to moderate 
the improvement at the longer L/D ratios. 

Both the absolute PCM temperature and the temperature dif- 
ference across the canister are seen to decrease as the number 
of elements is increased. Results shown in Figure 24 indicate 
that the LiF 22 CaF 2 , having a relatively low thermal conductivity, 
will require considerable sensible energy storage to ensure 
adequate potential to transfer heat to the heat pipe vapor during 
the eclipse phase. The increase in the number of canisters from 
15 to 30 serves to drop the peak temperature from about 12 0 above 
melting to about 3 0°. This is seen to occur in an exponentially 
decaying fashion. The minimum temperature is not shown here but 
will, by design, be quite close to melting. Thus the peak 
temperature above melting is essentially the same as the tem- 
perature difference. 



Figure 24. Effects of Element Number on Canister Temperature 

Difference. 

As seen in Figure 25 the effect of increasing the receiver 
length to diameter ratio is also to decrease both PCM temperature 
and temperature difference. In this case the temperature dif- 
ference for LiF 22 CaF 2 is decreased from about 100 to about 40 
degrees as the L/D ranges from 0.66 to 1.66. 
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Figure 25. Effects of Receiver L/D on Canister Temperature 

Difference. 



Figure 26. Effects of PCM Enrichment on Canister Temperature 

Difference. 

The PCM maximum and minimum temperatures are plotted against 
the element enrichment in Figure 26. As noted earlier, the effect 
of enrichment is to improve thermal conductivity but this occurs 
at the cost of increasing the length of the path through which 
heat conduction must occur. The utilization of 18% Boron Nitride 
or Graphite is to decrease the peak Element temperature from 
about 75 degrees above melting to around 50 and 30 respectively. 
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In this case the minimum temperature remains, by design, 
essentially at the melt temperature. The element temperature 
differences will, therefore, have the same value. 

3.3 Pebble Bed Receiver. 

The modeling of the pebble bed assumed that the pebble stacking 
would provide for efficient gas mixing so that the working fluid 
provides for a degree of PCM enhancement. This effect reduces 
the radial thermal resistance in the elements so that the impact 
of the maldistribution of input heat flux associated with low 
element circumferential pitch to diameter ratios is significantly 
reduced. As seen in Figure 27 the effect of varying the P/D 
between 1.4 and 2 is negligible on the minimum turbine inlet 
temperature, typically about 990 K throughout the range. The 
maximum turbine inlet temperature, corresponding to the end of 
solar incidence is seen to vary only between about 1050 down to 
1070 K over this same range. The effect of graphite enhancement 
is seen to be small but typically reduces the maximum temperature 
by 20 to 30 degrees over this same range. 



Figure 27. Effects of Element P/D on Turbine Inlet Tempera- 
ture . 

As seen in Figure 28, the effect of varying the receiver 
length to diameter ratio is also small. Varying its value from 
0.7 to 1.5 appears to decrease the turbine inlet temperature by 
about 15 degrees. This same variation occurs at both the peak 
and the minimum temperatures. 
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Figure 28. Effects of Receiver L/D on Turbine Inlet Tempera- 
ture. 



Figure 29. Effects of Element P/D on Canister Temperature 

Difference . 

The effect of element pitch to diameter ratio on the PCM 
temperature is shown in Figure 29. The effect is somewhat greater 
than that seen on the turbine inlet temperature. Here the P/D 
is varied from 1.4 to 2 the maximum PCM temperature at 54 minutes 
is seen to decrease from about 1122 to 1095 K. At the same time 
the minimum PCM temperature increases from about 1017 to 1028 K. 
The overall variation across the element is then seen to decrease 
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by about 25%. Both effects would be thought to be beneficial 
toward reducing high temperature, stress related structural 
problems . 

The PCM maximum temperature and temperature difference at 54 
minutes is also seen to be a weak function of the receiver length 
to diameter ratio as seen in Figure 30. As L/D varies from 0.75 
to 1.6 the peak PCM temperature is decreased from 1110 to 1100 
K. The temperature difference over this same range decreases 
from about 95 to 70 degrees. 



Figure 30. Effects of Receiver L/D on Canister Temperature 

Difference. 

3.4 Cavity Heat Pipe Receiver. 

As noted earlier each canister has associated with it a heat 
pipe wedge which serves effectively as a fin for thermal 
enhancement. Clearly the effect of additional fins will be to 
reduce temperature differentials within the element. Because 
the distances through which the energy is being driven is reduced 
with increased canisters, the temperature difference between the 
heat pipe fluid and the PCM melt temperature is reduced. The 
net effect is that cavity temperature variations between the end 
of the periods of solar incidence and solar eclipse are reduced. 
The requirement for small high temperature source variations is 
thought to be more severe for the Stirling engine than for the 
Brayton; it is noted that a Stirling receiver would therefore 
normally incorporate a higher number of elements than a corre- 
sponding Brayton receiver. 
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The results of the parametric study indicating the effect of 
number of canisters is shown in the Figure 31. Results indicate 
that temperature reductions average about 3.5 K per canister for 
the LiF 22 CaF 2 without enhancement and the BN enhanced PCM. The 
pyrolytic graphite demonstrates a somewhat, reduced effect with 
only about a 2.5 K reduction per canister. For the 100% PCM case 
this would imply that doubling the number of canisters from 18 
to 36 would decrease the cavity temperature variation from about 
175 to about 115 K. While this trend is somewhat reduced with 
the enhanced materials, 18% pyrolytic graphite will still undergo 
an improvement of from about 125 to about 65 K. The Ge, having 
an extremely high thermal diffusivity, shows a quite low tem- 
perature variation. In spite of the low value overall trends 
appear the same as for the remaining materials. Clearly the 
number of canisters is seen to be a primary factor in determining 
the magnitude of cavity temperature variation. 



Figure 31. Effects of Element Number on Canister Temperature 

Difference. 

A small apex angle will also tend to reduce temperature 
variations by reducing the maximum distance that heat must be 
transported into the solid. So long as a portion of the PCM 
remains partially melted the entire system is then coupled to 
the melting temperature. Short thermal distances are then seen 
to be close coupled whereas longer distances permit greater 
thermal swings. It is further noted that the small apex angles 
which tend to enhance heat transfer, involve configurations with 
large surface to volume ratios so that system weight will tend 
to suffer. Clearly additional trade-offs are required to optimize 
the system in this regard. 




The effect of apex angle on cavity temperature oscillations 
is shown in the Figure 32. The relationship appears to be nearly 
linear over the range shown, about 3 K per degree. Earlier 
results, not shown here, indicate that the enhancement effects 
begin to decrease above angles of about 22 degrees so that cavity 
temperature swings approach an asymptotic value. In effect 
decreasing the apex angle from the asymptotic limit to about 15 
degrees will only decrease temperature variations from about 125 
to a little over 100 K. A further reduction to a very sharp 
angle of 10 degrees will only drop the variation to 90 K. It is 
concluded that apex angle provides only modest improvement in 
cavity temperature variation and generally will be of secondary 
importance . 



Figure 32. Effects of Element Apex Angle on Cavity Tempera- 
ture Difference. 

Increasing element length will have the effect of decreasing 
the other element dimensions, specifically the depth and width. 
Since the latter two dimensions are the most significant in terms 
of the heat transfer problem, the net effect will be to decrease 
cavity temperature variations. 

As seen in Figure 33 for the unenhanced PCM, increasing the 
receiver length to diameter ratio from 0.5 to 2.0 produces a 
corresponding decrease in cavity temperature variations from 
about 160 to 80 K. For the two foam enhanced schemes the effect 
is similar; temperature variations are about halved. As would 
be expected the temperature variations for the metals, having 
very high thermal dif fusivities, are quite small. As with the 
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salts, the overall trend is toward decreased variations with 
increased L/D. In all cases the cavity temperature variation is 
seen to be highly dependent on receiver L/D. 



Figure 33. Effects of Receiver L/D on Cavity Temperature 

Difference. 

The effects of thermal enhancement are shown in Figure 34. 
The effect of enhancement tends to be nonlinear and small additions 
tend to produce limited results. This tendency is particularly 
seen for the BN case where additions of less than 10% show no 
significant overall improvement in thermal performance. A 20% 
addition has the overall effect of reducing cavity temperature 
swings from about 115 to about 70 K. The pyrolytic graphite 
would appear to offer more substantial thermal improvement; cavity 
temperatures are reduced to 90 and 50 K for foam additions of 10 
and 20%, respectively. Thermal enhancement is thereby seen to 
provide an effective means of further temperature control. 
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Figure 34. Effects of PCM Enhancement on Cavity Temperature 

Difference. 

3.5 Cavity Receiver. 

The variation in the turbine inlet temperature associated 
with the cavity receiver is shown in Figure 35. The range is 
seen to be relatively low with a minimum of 975 and a maximum of 
1015. This is the least variability of any of the designs 
considered. The effect is thought to be closely related to the 
large number of corrugations considered. Each of the other 
designs showed a tendency toward decreasin'*; variations with large 
number of elements. One advantage of this concept is that it 
seems to make the utilization of large numbers of PCM subdivisions 
without the weight and mechanical complexity of certain of the 
other designs. 

The effect of the large number of corrugations is to over- 
shadow the effects of other parameters on the turbine inlet 
temperature. Results showed little discernable variation in 
these temperatures with L/D, P/D enhancement, etc. Again, if 
the unit were designed with a number of corrugations more con- 
sistent with the number of elements in other receiver designs 
it is thought that this receiver would follow the established 
patterns. Because of the lack of variability in the parame- 
ters described the associated drawings have been omitted. 
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Figure 35. Turbine Inlet Temperatures Through Orbit. 

The effects of PCM enhancement on the PCM maximum and minimum 
temperatures are seen in Figure 36. Effects are seen to follow 
a pattern similar to that of other receivers. Clearly there is 
a tendency to reduce peak temperatures as enhancement is added; 
minimum temperatures appear to remain virtually unchanged. 
However the effect is seen to be relatively small. The receiver 
is thought to be designed with sufficient corrugations that the 
dimensions on individual PCM "elements" are small. The peak 
temperature is seen to decrease by about 20 degrees with 18% BN, 
about 35 degrees with a similar amount of graphite addition. 
Improvements in the temperature differentials are also relatively 
small, on the order of 20 to 30%. It should be recognized that 
with the transparency of the PCM much of the thermal energy 
bypasses the PCM reducing the required differential to produce 
melting in the PCM. 

A unique feature of this design is the ability to vary the 
energy distribution in the PCM by varying the absorptivity. This 
is done by controlling the quantity of absorbing material added. 
It would be anticipated that an optimal value of absorptivity 
could be found. A large value would tend to cause excessive 
heating at the PCM surface and result in increased temperature 
differentials. A low value would produce excessive absorption 
at the metal corrugation-PCM interface. Under an extreme condition 
virtually all energy would be added to the gas during periods of 
solar incidence with little storage occurring. Under these 
conditions no temperature difference would occur across the PCM. 
While the results shown in Figure 37 suggest that in the range 
indicated, the lower the absorptivity the better, this result is 
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differences increase with increasing L/D. This is an artifact 
of the method used in varying the L/D in the parameter study. 
The decision was made to hold the minimum receiver diameter 
constant while allowing the larger diameter to vary with length 
(The corrugation angles were held fixed) . The net result was to 
increase the PCM radial dimensions as length was added. 

It may be observed that the minimum PCM temperature is held 
slightly below the melting point at all times so that this 
parameter shows little variability. The maximum temperature 
increases relatively rapidly as the radial PCM dimension 
increases. (Note that the added fraction of absorbing material 
must be continually adjusted to maintain a constant absorptivity. ) 
Under these conditions the temperature differential is seen to 
nearly double, from 80 to 145 degrees as the L/D is increased 
from 0.66 to 1.33. 



Figure 38. Effects of Receiver L/D on Canister Temperature 

Difference. 

As seen in Figure 39 the effect of increasing the number of 
corrugations is similar to that of increasing the number of 
elements for other concepts. Here the range extends from 75 to 
150. Under these conditions the minimum temperature shows little 
variation, remaining slightly below the melting point. The 
maximum temperature decreases from about 1165 to about 1130 for 
the unenhanced PCM. Additions of 18% BN decreased this peak 
temperature from 1135 to 1110 over this same range. Graphite 
provides for a further reduction ranging from 1110 to 1090 K. 
From the increasing slope of the curve it would appear that 











4 Discussion 

It should be emphasized that in a evaluation of this type, 
including only thermal performance of the relative concepts, relative 
comparisons are somewhat limited and should not be used as a sole 
basis for receiver selection. Moreover there remain a certain 
arbitrariness in establishing certain of the configurations. 
Specifically those designs including canisters were not credited 
with the thermal enhancement properly attributed to that design. 
This was done in order to illustrate a continuity in the effects 
of thermal enhancement where it was assumed that PCM distribution 
would be accommodated by the foam filler itself. Further it is not 
clear to the authors that the effects of potential liquid PCM 
circulation and thermal insulation due to possible void coalescence 
and migration have been equitably evaluated between competing 
concepts. Finally the use of the CAV2 generated incident heat 
fluxes may unfairly penalize the non heat pipe concepts by failing 
to account for radiant redistribution of the incident heat flux 
within the receiver. Nevertheless it is believed that the results 
presented are useful in a number of ways. In a general sense they 
may serve to provide reasonable estimates of the performance 
achievable. They also provide an indication of the benefits 
attributable to various enhancement alternatives. Finally they 
serve to indicate that range of individual parameters which tend 
to optimize thermal performance within each of the respective 
designs . 

4 . 1 Attainable Performance 

Attainable thermal performance is evaluated here in two 
distinct ways. The first involves the variation in receiver 
outlet temperature, generally termed either turbine inlet tem- 
perature, heat pipe temperature, or engine head temperature. The 
second involves PCM peak temperatures or PCM thermal differences. 
These parameters relate to the thermal stress and thermal cycling 
within the PCM containment and have a direct bearing on the 
ultimate reliability of the receiver. 

The cyclic thermal variation for the base line receiver is 
shown in Figure 12 where the overall variation is seen to be of 
the order of 60 K. The basic problem of optimizing performance 
with this design is the maldistribution cf incident energy in 
the axial direction. Excess incident energy is provided at the 
aperture end of the receiver tending to provide excessive sensible 
heating of the PCM during the latter stages of the heating cycle. 
Insufficient heating occurs at the opposite end of the receiver 
so that complete solidification and sensible heat removal during 
the eclipse phase remains a problem. The design problem then 
becomes one of making a trade off between excess temperature near 



the receiver inlet versus increased engine hot side temperature 
variations. The most direct solution of the problem is to increase 
the amount of PCM in the receiver. 

The corresponding plot for the heat pipe receiver is shown 
in Figure 20. It is important to note here that the number of 
elements in the heat pipe system is reduced to 20 from the original 
75 used in the base line receiver. It is also important to note 
that the calculated results for the heat pipe system do not 
include any enhancement effects from the canister side walls. 
As a consequence the LiF 22 CaF 2 results indicate somewhat larger 
temperature swings than would be observed in such a system. Based 
on prior runs a reasonable approximation for such a system would 
be to assume that the effects of adding 12% metal fins would be 
to produce about 2/3 of the improvement of adding 18% BN. Thus 
the anticipated temperature variation for the system is thought 
to be about 94 K. Increasing the number of elements from 20 to 
30 with the heat pipe system, as seen in Figure 22, would serve 
to decrease the temperature difference to about 50 K. The two 
systems are strikingly similar in terms of individual components 
and would seem to offer the most viable direct comparison of any 
two receivers under consideration. The improvements associated 
with the heat pipe configuration may be attributed to a combination 
of factors. First the heat pipe serves to smooth the incident 
heat flux. Second, energy utilized by the engine during the 
period of solar incidence is not forced through the PCM, but 
rather is transported directly from the aperture section of the 
heat pipe directly to the engine. The resultant reduction in 
heat flux through the PCM can only serve to reduce the required 
temperature differential. 

The engine hot side temperature variation for the cavity 
receiver is shown in Figure 35. This design also utilized a 
direct solar incidence on the PCM with the working fluid flowing 
past the PCM. As with the base case this design has the problem 
of excessive heat addition near the aperture and insufficient 
heat flux at the outlet end. It does benefit from a somewhat 
higher number of canisters than the base line receiver and from 
a truncated cone shape which redistributes the incident heat flux 
to an extent. The result of these two improvements is seen in 
reducing the working fluid temperature variation to about 47 K. 

It is concluded that an overall temperature variation of 
between 50 and 100 K is reasonably achievable with an unenhanced 
LiF 22 CaF 2 PCM receiver. As will be noted in the following 
sections, the lower range of values will typically be associated 
with arrangements optimized for thermal performance rather than 
weight . 
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The range of radial temperature variation across the indi- 
vidual elements is also of importance, not to engine performance, 
but rather as an indication of the thermal stresses established 
within the PCM containment. This variation is normally a maximum 
at the end of solar incidence and is therefore evaluated at that 
time. Shown in Figure 19 is the maximum and minimum temperatures 
within that range of canisters nearest the aperture of the base 
line case. Since the reference condition is at a P/D of about 
1.66 it can be seen that the maximum and minimum temperatures 
correspond to about 172 and 148 degrees above melting or 1211 
and 1187 respectively. Note that the variation across the element 
is only 24 degrees but the absolute values are relatively large. 
This is due to the high axial incident flux near the aperture. 

The cavity receiver, as shown in Figure 39, tends to have 
slightly improved PCM temperature differentials over the base 
line case. With 125 corrugations the PCM maximum is 1132 and 
minimum is 1022 K. The reduction in the maximum temperature can 
be attributed partially to the improved incident flux distribution 
associated with the truncated cone configuration. It is also 
important to note that the coolant passages in this design are 
suitable for conforming to axial needs. In this analysis it has 
been assumed that the convection coefficient is axially adjusted 
to provide an axial NTU distribution which matches the fraction 
of heat removal to the fraction of incident energy. Upon reviewing 
the final results it was realized that this effect could be 
carried to a greater extent. That is the matching of the rate 
of heat removal to the incident heat flux leads to a situation 
in which the orbital temperature variation tends to produce axial 
PCM temperature patterns in which the minimum temperature appears 
as something of a mirror image of the maximum temperature with 
an axis about the melting point. A further refinement of the 
method would increase the rate of heat removal at axial locations 
corresponding to the highest PCM temperatures. This would distort 
the mirror image decreasing the peak temperatures and producing 
significant sensible heat removal during the heat removal phase 
at these same locations. 

In comparison both the heat pipe and the cavity heat pipe 
receivers operate at maximum and minimum PCM temperatures of 
about 60 and 0 degrees above melting or 1099 and 1039 K, 
respectively. The magnitude of the temperatures is significantly 
reduced due to smoothing of the axial flux in the heat pipe. The 
temperature differential is, however, considerably larger. This 
may be attributed to the reduced number of elements and the 
subsequent increase in individual canister radius. 
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The pebble bed arrangement is by it nature somewhat different 
in terms of element temperature variations. The values shown in 
Figure 30 indicate that the temperature difference across the 
PCM mass ranges from 1025 to 1110 K. Note that this difference 
is not across a single sphere but represents the temperature 
difference between spheres at the center and edge of the bed. 
For this reason comparisons of metal stress between these elements 
and the alternative designs are difficult. In this design the 
PCM sphere size is essentially 1/3 of the element containment 
ID. A simple analysis would then suggest a temperature difference 
across each sphere of about 30 to 35 degrees. 

4 . 2 PCM Thermal Enhancement 

The effects of thermal enhancement on the cycle hot side 
temperature are shown for the base line receiver in Figure 13. 
For the heat pipe receiver these same results are shown in Figure 
20. Neither case provides a clear comparison of the various 
enhancement schemes. As noted earlier, the set of conditions 
chosen for the base liner receiver do not lend themselves to the 
use of Ge or 18% additions of BN or pyrolytic graphite. As a 
consequence consideration has been limited to 10% additions of 
the latter two enhancement materials. Moreover both cases are 
somewhat complicated by the question of whether and how side wall 
should be considered. The approach proposed here is simply to 
compare the results neglecting the side walls but to weight 
certain of the other cases more heavily in making the final 
comparison. 

In the base line case the cyclic variations are reduced from 
about 77 K to 43 and 37 K for a 10 % addition of BN and graphite, 
respectively. This corresponds to a 44% and a 52% reduction in 
variation from the base line. The enhancement of the heat pipe 
receiver is seen to reduce the variation from about 110 K to 
about 75, 60 and 5 K for the 18% BN 18% graphite and the Ge, 
respectively. This results in a 32% improvement for the BN, and 
a 45% improvement for the graphite. Note that while the while 
the enhancement fraction is much lower for the base line case 
the improvement achieved is greater. The introduction of Ge as 
a phase change material reduces the overall temperature variation 
to near zero and would clearly be the most advantageous choice 
from thermal considerations. 

A somewhat simpler comparison is made in the various 
enhancement schemes in the cavity heat pipe receiver. Referring 
to Figure 31, it is seen that using the reference 36 elements 
the introduction of thermal enhancement reduces the variation 
from 115 to 90 and 62 K, respectively. The improvement is seen 
to be a reduction of 22% and 45% respectively. While the Ge PCM 
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has not been evaluated for this same number of elements, certainly 
the variation will be less than for the 2* element case or less 
than 15 K. The corresponding reduction is then greater than 87%. 

The cavity receiver provides a second relatively direct 
comparison. As seen in Figure 34, the BN serves to reduce the 
variation from about 117 to 95 K, or about 19%. Similarly the 
graphite reduces the variation to about 70 K or about 40%. 

Summarizing the cases above it is concluded that thermal 
enhancement with BN would normally be expected to decrease engine 
hot side temperature variations to between 20 and 25%, graphite 
to between 40 and 45% and Ge to between 85 and 95%. 

Canister temperature differences are shown for the cavity 
heat pipe receiver can be taken from Figure 32 recognizing that 
at 54 minutes the minimum temperature will be at the back wall 
of the element where only a residual solid mass of PCM remains. 
As a consequence the element temperature difference will be 
between the melting temperature and the peak temperature. It 
would appear then that the addition of 18% BN would serve to 
reduce temperature variation across the element, from 65 to 48 K 
or by 26%. Graphite enhancement reduces the differential to 32 
K or by 51%. Germanium PCM would encounter only about a 5 K 
difference at this same time representing a reduction of 92%. 

The cavity receiver indicates a much larger temperature 
difference across the elements as indicated in Figure 36. Here 
the unenhanced material operates with a peak temperature dif- 
ference of about 120 K. Eighteen percent additions of BN reduce 
this to about 100 K or about 17% and similar additions of graphite 
to about 79 K or by 34%. 

Summarizing these conclusions it appears that BN enhancement 
would potentially reduce peak PCM temperature differences by 
between 15 and 30%. Similar additions of graphite would reduce 
these same differences by between 35 and 50%. Germanium PCM 
would virtually eliminate such differences, reducing them by 
about 90% as compared to the unenhanced salt. 

4.3 Parametric Studies 

In reviewing the thermal performance of each of the design 
concepts, it becomes clear that most of the corresponding per- 
formance curves appear similar in general form. That is the 
cycle hot side temperature, whether it be termed the turbine 
inlet temperature or heat pipe temperature, tends to become more 
uniform as the effective thickness of the PCM is decreased. 
Several methods remain for controlling the effective thickness 
in each of these designs. 



Generally, increasing the receiver length will tend to 
redistribute PCM so as to provide a smaller radial path for energy 
storage. In the study shown here this effect is shown as the 
receiver length to diameter (L/D) ratio. An apparent exception 
to decreased temperature variations with increased L/D is the 
case of the cavity receiver. This case involves the use of an 
absorbing phase in the PCM. It is thought that the higher 
temperatures are associated with the steeper energy gradient 
associated with the increased concentration of absorbing material 
utilized to maintain a constant average absorptivity. 

In each of the designs it was found that increasing the number 
of elements, or corrugations in the case of the cavity receiver, 
tended to reduce variations in the cycle hot side temperature. 
This again is attributable to the reduction in the effective 
thickness of PCM through which the heat flux must travel. In 
those cases which did not utilize the heat pipe concept to 
redistribute the incident heat flux, it was found that increasing 
the element pitch to diameter ratio would tend to redistribute 
incident heat flux in a more uniform method also reducing tem- 
perature differences. 



5 Conclusions. 

While the base line concept represents a relatively high degree 
of development in thermal receiver design, substantial improvements 
in thermal performance may still be achievea. Generally these are 
accomplished through ' one of several methods including (1) the use 
of heat pipes to redistribute incident heat flux in a uniform manner, 
(2) better utilization of the PCM material with less sensible heating 
and cooling, (3) introduction of concepts which completely or 
partially bypass energy flow around the PCM and directly to the 
engine, (4) introduction of concepts which can practically accom- 
modate larger number of elements or introduce other effective methods 
of reducing the thickness of PCM through which heat transport occurs 
or (5) introduction of enhanced thermal conductivity materials. 

Within this range of design features it is possible to reduce 
peak PCM temperature from about 1210 to around 1090 K without the 
benefit of thermal enhancement. The base line receiver is found 
to provide relatively uniform engine hot side temperatures with 
variations of about 60 K. While certain of the concepts may in 
fact increase this value, reductions to about 45 K are possible 
while retaining the LiF 22 CaF 2 ' PCM. 

The introduction of BN enhancement is found to marginally 
beneficial in terms of improved thermal performance. The development 
of a graphite enhancement material would provide somewhat higher 
levels of performance and would receive a much higher priority. 
The development of Germanium or other metallic PCM offers a quite 
high level of thermal performance is, if containment or other 
potential development problems can be overcome, offer substantial 
improvements over all other enhancement schemes. 
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APPENDIX A 


PROGRAM FOR THE THERMAL ANALYSIS 
OF THE 

BASE LINE RECEIVER 
AND THE 

HEAT PIPE RECEIVER 



A1 





C PROGRAM SOLDTESS 

C 

C SOLAR DYNAMIC TRANSIENT ENERGY STORAGE SIMULATION 

C 

C This program is the main routine for the TESS code. 

C Its purpose is to read & print data and to call 
C calculational subroutines so as to complete the thermal 

C analysis. Checks for convergence are made here. 

C 

c************************************************************- 

c 

C************************************************************' 

c 

C NOMENCLATURE 


C 

C System dimensions 

C 

C CS = Coefficient for Simpson 1/3 Integration 

C DR = Diff. element dimension in radial direction, m 

C DPHI = Diff. element dim. in azimuthal direction, rad 

C DMAX = Maximum move in melt position in one iteration 

C DMELT = Fractional movement of melt position in element 

C DTIME = Time increment, sec 

C DZ = Diff. element dimension in axial directrion, m 

C GASMAS = Coolant flow rate, kg/sec 

C LENGTH = Axial length of stack, m 

C OVRRLX = Over relaxiation factor on melt movement 

C POS = Axial position of current node, m 

C R = Radius of element (I,J), m 

C RI = Outside radius of inner can ring, m 

C RICAN = Inside radius of inner can ring, m 

C RMELT = Fraction of PCM in element melted 

C RN = Radius at outside edge of element (I,J), m 

C RO = Inside radius of outer can ring, m 

C ROCAN = Outside radius of outer can ring, m 

C RS = Radius at inside edge of node (I,J), m 

C ROMELT = Previous fraction of PCM melted 

C T = Projected relative temperature of node 

C TCI J = Absolute coolant temperature of gas at pos., K 

C TGAS = Relative coolant temperature at axial position 

C TGAS IN = Absolute coolant temp, entering axial position, K 

C TGASP = Absolute coolant temperature at mid DZ , K 

C THGAS = Absolute coolant temperature at inlet, K 

C THETA = Absolute temperature of node (I,J), K 

C THETAG = Absolute temperature of gas, K 

C THMELT = Absolute melting temperature of storage media, K 



A2 



TIME = Time into orbit run, min 
TIMER = Time of orbit, min 

TO = Previous relative temperature of node (I,J) 

TP = Projected relative temperature of node (I,J) 

WALLS = Ratio of wall conduction area to storage area 

Material properties 

CPGAS = Coolant specific heat, kJ/kg K 
CVLIQ = Liquid specific heat, kJ/kg K 
CVMET = Ring specific heat, kJ/kg K 

CVPURE = Specific heat of PCM (solid & liquid) kJ/kg K 
CVSOL = Solid specific heat, kJ/kg V 
DENLIQ = Liquid density, kg/m**3 
DENMET = Ring density, kg/m* *3 

DENPURE = Density of PCM (solid & liquid) kg/m**3 
DENSOL = Solid density, kg/m**3 
F = Smooth tube adiabatic friction factor 
GASMIX = Fraction of gas 1 in coolant mix 
H = Convection coefficient of gas W/m**2 K 
HFL = Heat of fusion, kJ/kg 

K = Eff. thermal conductivity of node (I,J), W/m K 
KE = Thermal conductivity of node (I,J+1), W/m K 

KN = Thermal conductivity of node (I+1,J), W/m K 

KP = Thermal conductivity solid with void, W/m K 
KS = Thermal conductivity of node (I-1,J), W/m K 

KW = Thermal conductivity of node (I,J-1), W/m K 

K_1 = Thermal conductivity constant (See GASPROP) 

K_2 = Thermal conductivity constant (See GASPROP) 

K_3 = Thermal conductivity constant (See GASPROP) 

K_4 = Thermal conductivity constant (See GASPROP) 

KLIQ = Liquid thermal conductivity, W/m K 
KMET = Ring thermal conductivity, W/m K 

KPURE = Eff. thermal cond. of PCM (solid & liq.) W/m K 
KSOL = Solid thermal conductivity, W/m K 
MW_ = Molecular weight of gas component _ 

VISC_1 = Viscosity constant (see GASPROP) 

VISC_2 = Viscosity constant (see GASPROP) 

VISC_3 = Viscosity constant (see GASPROP) 

VISC_4 = Viscosity constant (see GASPROP) 

VOID = Void fraction in element 

Indicies 

I = Counter in radial direction 
ICOUNT = Loop counter 

I FLAG = Flag to indicate insufficient convergence 

IMAX = Maximum number of loops 
J =. Counter in azimuthal direction 
L = Counter in axial direction 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c* 

c 

c 


c 


LIO = Counter on output file 

NORBIT = Number of orbit cycles 
NI = Number of radial nodes 
NJ = Number of azimuthal nodes 
NL = Number of axial nodes 

NTIME = Number of time increments per orbit 
Matrix coefficients 

A = Conduction matrix coefficient on 1-1 node 
B = Conduction matrix coefficient on I node 
C = Conduction matrix coefficient on 1+1 node 
D = Conduction matrix coefficient on constant vector 

Thermal resistances & capacitances 

AE = Thermal conductance of J-l node, W/K 

AN = Thermal conductance of I+l node,. W/K 

AS = Thermal conductance of 1-1 node, W/K 

AW = Thermal conductance of J+l node, W/K 

APO = Thermal capacitance of node (I, J) , kJ/m K 

Heat fluxes 

QAZ = Outside surface heat flux, W/m**2 
QDIR = Direct solar heat flux, W/m**2 

QGASP = Total heat added to gas at position DZ/2, kJ 
QGASDZ = Total heat added to gas at postion DZ , kJ 
QIJ = Radiation to element (I,J), W/m**2 
QIDIR = Indirect solar heat flux, W/m**2 
QINTR = Heat input on internal surface, W/m**2 
QZ = Integrated avg. axial multiplier on heat flux 
QZZ = Local axial multiplier on heat flux 
QZ_ = Constant on fit of QZ 

* 

***************************************************************** 


PROGRAM SOLDTESS 

REAL KLIQ,KMET,KSOL, LENGTH, MIN (11) , NTU 
DIMENSION FLOGAS (10) , QDIR (10) ,QINDIR(10) ,QINT(10) 

DIMENSION TGAS (10) , THETA (10) , THGAS ( 10 , 10 , 10 ) , TP ( 15 , 10 , 10 ) 

COMMON/ BLOKA/CVLIQ , CVMET, CVSOL, DENLIQ, DENMET, DENSOL, HFL, 

, KLIQ , KMET , KSOL , LENGTH , PI , RI , RICAN , RO , ROCAN , 

,RMELT(15, 10, 10) , R0MELT ( 15 , 10 , 10) , VOID, WALLS 
COMMON/ BLOKB/ F, GASMIX 

COMMON/ BLOKC/ DR, DTI ME , DPHI , DZ , HEAT, QAZ , QZ , TO ( 15 , 10 , 10 ) 
COMMON/ BLOKD/CPG AS , DELTDZ , DELTP , GASMAS , H , NTU , TGAS IN , TGASP , 
, THMELT 
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COMMON/ BLOKE/A ( 15) , B ( 15) , C ( 15) , D ( 15) , T ( 15 , 10 , 10) 

COMMON/ BLOKF/ I, J, L,NI,NJ,NL 
COMMON/ BLOKG/POS , QZ1 , QZ2 , QZ3 , QZ4 
C 

C*************************************************************** 

c * 

C OPEN INPUT AND OUTPUT FILES * 

C * 

C*************************************************************** 

c 

OPEN ( 8 , FILE= ' INPUT . TXT ' ) 

OPEN ( 7 , FILE= ' OUTPUT . TXT ' ) 

11=1 

PI=3 . 14159 

READ (8,2100, ERR=9 3 ) IMAX , NORBIT , NI , NJ, NL, NT I ME 
READ (8,2200, ERR=9 4 ) RICAN , RI , RO , ROCAN , LENGTH , DUM1 
READ (8,2200, ERR=9 5 ) CVLIQ , DENLIQ , KLIQ , WALLS , THMELT , DUM2 
READ (8,2200, ERR=9 6 ) CVSOL , DENSOL , KSOL , ERRMAX , HFL, DUM3 
READ (8,2200, ERR=9 7 ) CVMET , DENMET , KMET , GASMI X , TIMER , BORL 
READ (8,2200, ERR=98 ) QZ 1 , QZ2 , QZ3 , QZ4 , SORL 
1 READ ( 8 , 2200 , ERR=99 , END=5 ) MIN ( II ) ,TGAS(II) ,FLOGAS(II) , 

, QDIR ( II ) ,QINDIR(II) ,QINT(II) 

11=11+1 
MIN ( II ) =1 . E9 
GO TO 1 
C 

C * 

C INITIALIZE LOOP COUNTERS AND ERROR FLAGS * 

C * 

C*************************************************************** 

c 

5 DELTP=0 . 0 
ITEMP=0 

DR= (RO-RI ) /FLOAT (NI-2 ) 

DPHI=PI/ FLOAT (NJ) 

DZ=LENGTH/ FLOAT ( NL) 

DTIME=60 . 0*TIMER/FLOAT (NTIME) 

C 

DO 10 11=1, NI 
DO 10 JJ=1 , NJ 
DO 10 LL=1 , NL 

T ( II , JJ , LL) =-50 . 0 
TP (II, JJ,LL)=T(II, JJ,LL) 

C 

IF (II .LT. NI ) THEN 

RMELT ( II , JJ , LL) =0 . 0 
ELSE 

RMELT ( II , JJ , LL) =1 . 0 
ENDIF 



non 
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C 

10 CONTINUE 

C' *************************************************************** 

* 

BEGIN CALCULATIONAL LOOPS ON LENGTH, TIME, ANGLE * 

* 

C' *************************************************************** 
C 

DO 100 NO=l , NORBIT 
NPUT=0 
TIME=0 . 0 
C 

DO 100 N=1 , NTIME 

IF (MIN (NPUT+1) .LE. TIME+0 . IE-6 ) NPUT=NPUT+1 
GASMAS=FLOGAS (NPUT) 

THGAS ( N , NO , 1 ) =TGAS ( NPUT ) 

TIME= FLOAT (N) *DTIME/60 . 0 
C 

DO 100 L=1 , NL 

POS=DZ * FLOAT (L) 

TGAS IN=THGAS ( N , NO , L ) 

C 

DO 20 11=1, NI 
DO 20 JJ=1 , NJ 

R0MELT ( II , JJ , L) =RMELT ( II , JJ , L) 

TO ( II , JJ , L) =T ( II , JJ , L) 

20 CONTINUE 

C 

ICOUNT=0 
30 IFLAG=0 

ICOUNT=ICOUNT+l 
TGAS P=THGAS ( N , NO , L ) + DE LT P 
C 

DO 40 J=1,NJ 

CALL GASPROP (RICAN, T ( 1, J,L) , DZ) 

CALL SOLAR ( LENGTH, QDIR( NPUT) , QINDIR (NPUT) ,ROCAN) 
CALL CONDUCT (BORL, SORL, QINT (NPUT) ) 

CALL SOLVER 

CALL HEATGAS ( T ( 1 , J , L) ) 

40 CONTINUE 
C 

DO 50 11=1, NI 
DO 50 JJ=1,NJ 

ERROR=ABS (T(II,JJ,L) -TP(II,JJ,L) ) /THMELT 
IF (ERROR .GT. ERRMAX) IFLAG=2 
IF (ERROR .GT. ERRMAX) ITEMP=ITEMP+1 
TP (II , JJ , L) =T (II , JJ, L) 

50 CONTINUE 

C 

IF ( I COUNT .GT. IMAX) THEN 




non 


WRITE (7, 1000) ITEMP 
IFLAG=0 
ENDIF 
C 

IF ( IFLAG .GE. l)GO TO 30 
THGAS ( N , NO , L+ 1 ) =THGAS (N,NO, L) +DELTDZ 
IF (NO .LT. NORBIT) GO TO 100 
C 

WRITE (7, 1050) 

WRITE(7, 1250) 

WRITE (7,1100) POS , TIME , NO 
WRITE (7, 1150) QZ, HEAT 
DO 70 1=1, NI 

DO 60 J=1,NJ 

THETA ( J) =T ( I , J , L) +THMELT 
60 CONTINUE 

WRITE(7, 1200) (THETA (J) ,J=1,NJ) 

70 CONTINUE 

C 

WRITE (7, 1300) TGASP, THGAS (N, NO, L+l) ,H 
C 

WRITE (7, 1050) 

WRITE (7, 1250) 

DO 80 1=1 , NI-1 

WRITE (7, 1350) (RMELT ( I , J , L) , J=1,NJ) 

80 CONTINUE 

WRITE (7, 1350) (RMELT ( 1 , J , L) , J=1,NJ) 

C 

100 CONTINUE 
C 

C* ***************************************** * ******************* 

* 

INPUT AND OUTPUT FORMAT STATEMENTS * 

* 

C* ************************************************************** 
C * 

1000 FORMAT (IX, 'FAILED TO CONVERGE AFTER ',16,' ITERATIONS.') 
1050 FORMAT (IX,' ') 

1100 FORMAT ( IX, 'Axial Position = ',1F5.3,' m. Time = ', 

& 1F6.1,' min. Cycle = ',11,'.') 

1150 FORMAT ( IX, 'Heat Input Frac. = ',1F6.4,' Heat = ',1F9.4) 
1200 FORMAT ( 10F7. 1) 

1300 FORMAT ( IX, 'Avg. Gas Temp. = ',F6.1,'K Out. Gas temp.', 

&F6 . 1 , ' K H = ' ,F7.4, 'KW/m2-K') 

1350 FORMAT (10F7. 4) 

1400 FORMAT (IX, 'INPUT FORMAT ERROR IN 1ST LINE') 

1401 FORMAT (IX, 'INPUT FORMAT ERROR IN 2ND LINE') 

1402 FORMAT (IX, 'INPUT FORMAT ERROR IN 3RD LINE') 


1403 

1404 

1405 

1406 
2100 
2200 

93 

94 

95 

96 

97 

98 

99 


FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

STOP 

WRITE ( 

STOP 

WRITE ( 

STOP 

WRITE ( 

STOP 

WRITE ( 

STOP 

WRITE ( 

STOP 

WRITE ( 

STOP 

WRITE ( 

STOP 

END 


(IX, 'INPUT FORMAT ERROR IN 4TH LINE') 

(IX, 'INPUT FORMAT ERROR IN 5TH LINE') 

(IX, 'INPUT FORMAT ERROR IN 6TH LINE') 

(IX, 'INPUT FORMAT ERROR IN TRANSIENT SECTION') 

/ T r r i a atvi a a \ 


(I6,5I10,2F10.0) 
( 2F10 . 4 , 4E10 . 4 ) 

*,1400) 

*,1401) 

*,1402) 

*,1403) 

*,1404) 

*,1405) 

*, 1406 ) 
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C************************************************************** 


c * 

C SUBROUTINE SOLAR * 
C * 
C This program calculates the solar incidence on * 
C the outside perimeter of the containment ring * 
C * 
C q" = q" (direct) *cos (angle) + q" (indirect) * 
C * 


c************************************************************** 

c************************************************************** 

c 

SUBROUTINE SOLAR ( LENGTH , QDIR , QINDIR , ROCAN) 

C 

DIMENSION CS (21) 

REAL LENGTH 

COMMON/ BLOKC/ DR , DTIME , DPHI , DZ , HEAT , QAZ , QZ , TO ( 15 , 10 , 10 ) 
COMMON/ BLOKF/ I , J , L , NI , N J , NL 
COMMON/ BLO KG/ POS , QZ1 , QZ2 , QZ3 , QZ4 

DATA CS(1) , CS ( 3 ) ,CS(5) ,CS(7) ,CS(9) ,CS(11) ,CS(13)/ 

/l. 0,6*2. 0/ 

DATA CS ( 15 ) , CS ( 17 ) , CS ( 19 ) , CS ( 2 ) , CS ( 4 ) , CS ( 6) , CS ( 8 ) / 

/3*2.0, 4*4.0/ 

DATA CS(10) , CS ( 12 ) ,CS(14) ,CS(16) ,CS(18) ,CS(20) ,CS(21)/ 
/6*4. ,1./ 

DATA PI/3.14159/ 

C 

C* ************************ * * * * * * * * * * *************** * * * * * * * * * * * * * 


c * 

C QZ USES SIMPSON'S 1/3 INTEGRATION TO FIND THE * 
C INTETRATED AVERAGE VALUE OF THE AXIAL POWER CURVE * 
C OVER THE INTERVAL DZ * 
C * 


C*************************************************************** 

c 

XL= ( POS-DZ ) /LENGTH 
QZ = 0 . 0 
C 

DO 15 11=1,21 

IF (II .GT. 1) XL=XL+0. 05*DZ/LENGTH 
QZZ=QZ 1*XL*EXP (1.0- (XL/QZ2 ) **QZ3 ) +QZ4 * (XL) **0.5 
QZ=QZ+CS (II)*(QZZ)/60.0 
15 CONTINUE 
C 

C*********************************************-*-**************** 

c * 

c FIND THE PROJECTION OF THE DIRECT SOLAR INCIDENCE * 

C * 

C ***** ****************************** *************************** 
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c 

ANGLE1=PI*FL0AT(J) /FLOAT (NJ) 

ANGLE2=PI* FLOAT (J-l) /FLOAT (NJ) 

C 

IF ( J+J-2 .GT. NJ) THEN 
' SINANGLE=0 . 0 
ELSEIF (J+J .GT. NJ) THEN 

SINANGLE=1 . O-SIN (ANGLE2 ) 

ELSE 

S INANGLE=S IN ( ANGLE 1) -SIN (ANGLE2 ) 

ENDIF 

************************************************************** 

* 

ADD DIRECT AND INDIRECT INCIDENT SOLAR RADIATION * 

* 

C********************** ****************** ********************** 
C 

QAZ=(QDIR*SINANGLE + QINDIR*DPHI) *ROCAN*DZ*QZ 
IF ( J .EQ. 1 ) HEAT=0 . 0 
HEAT=HEAT+2 . 0*QAZ 
C 

RETURN 

END 
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Q***** ****************************** ******* k ******************* 
Q* ************** * ********************************* ************* 

C * 

C SUBROUTINE GAS PROP * 

C * 

C This program calculates the effective thermal * 

C properties and convective coefficients for the gas coolant * 
C The gas coolant is assumed to be a mixture of He and Xe. * 
C Properties for these two constituents are included as * 

C curve fits with the corresponding coefficients included * 
C as data within the subroutine. * 

C * 

Qk k k * -k k * *k k * * * * k * * * * * -k k k * * k * * * * * * * * * * * * * * * * * k A * * * * * * * * * * * * k * k * * ★ 
Qk k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k 


c 

Qkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

c * 

C REFERENCES * 

C M. R. Vanco, "Analytical Comparison of Relative Heat * 

C Transfer Coefficients and Pressure Drops of Inert * 

C Gases and Their Binary Mixtures", Lewis Research Center* 

C * 

C A. J. Chapman "Fundamentals of Heat Transfer" * 

C Macmillan, 1987, p 329. * 

C * 

Qk k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k 


Qk k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k 

c * 

C EQUATIONS * 

C * 

C K1 = Kll + K12 * (TGASP - K13)**K14 * 

C K2 = K21 + K22* (TGASP - K23)**K24 * 

C * 

C VISC1 = VISC11 + VISC12* (TGASP - VISC13 ) **VISC14 * 

C VISC2 = VISC21 + VISC22* (TGASP - VISC23 ) **VISC24 * 

C * 

Q************************************************************** 

c 


SUBROUTINE GASPROP (RICAN , T1J , DZ ) 

C 

REAL K, K1 , K2 , Kll , K12 , K13 , K14 , K21 , K22 , K23 , K24 , KGAS , 

+MW , MW1 , MW 2 , NTU , NU 
C 

COMMON/ BLOKB/ F , GASMIX 

COMMON/ BLOKD/CPGAS , DELTDZ , DELTP , GASMAS , H , NTU , TGAS IN , TGASP , 
, THMELT 
C 

DATA GAMMA, PI , RUNIV, MW1 , MW2/1 .66,3 . 14159,8 . 3 14 3 4 , 4 . , 13 1 . / 
DATA Kll, K12,K13,K14/0. 1580E-3,0. 9086E-6, 333 .333, 0.794/ 





non 


All 






DATA K21 , K22 , K23 , K24/0 . 6404E-5 , 0. 04 54 E-6, 333. 333, 0.800/ 
DATA VISC 11 , VISC12 , VISC13 , VISC14/ 

/ . 2 158E-4 , . 1022E-6, 333 . 3 , .825/ 

DATA VISC2 1 , VISC22 , VISC2 3 , VISC24/ 

/ . 253E-4 , . 2452E-6, 333.333, .77/ 

C 

FI (X, Y) = (1.0+2.41* (X-Y) * (X-0 . 142*Y) / (X+Y) **2 ) 

F2 (XI , X2 , X3 , X4 , X5 , X6 ) = 

=( (1.0+ (X1/X2) **0.5* (X3/X4) **0.25) **2)/ 

/(2. 0*2. 0**0. 5*(1. 0+X5/X6) **0.5) 

C 

C*********************************************************** 

* 

FIND THE THERMAL CONDUCTIVITY & VISCOSITY OF PURE GASES * 

★ 

0 *********************************************************** 

c 

K1=K11+K12* (TGASP-K13 ) **K14 
K2=K21+K22* (TGASP-K23 ) **K24 

VISC1=VISC11+VISC12 * (TGASP-VISC13 ) **VISC14 
VISC2=VISC21+VISC22* (TGASP-VISC23 ) **VISC24 
VISCW1=VISC11+VISC12* (T1J+THMELT-VISC13) **VISC14 
VISCW2=VISC2 1+VISC22* (T1J+THMELT-VISC23 ) **VISC24 
C 

C **** ********************************************* ********** 

C * 

C FIND THE WEIGHTING PARAMETERS FROM REFERENCE * 

C * 

C*********************************************************** 

C 

PHI 12=F2 ( VISC1 , VISC2 , MW 2 , MW1 , MW1 , MW 2 ) 

PHI2 1=PHI12* (VISC2*MW1) / (VISC1*MW2) 

PSI 12=F2 ( K1 , K2 , MW1 , MW2 , MW1 , MW2 ) *F1(MW1,MW2) 

PSI21=F2 (K2 , K1 , MW2 , MW1 , MW2 , MW1) *F1(MW2,MW1) 

C 

Qlc it ic Jc * ic ic ic 1c ic ic 1c ic ie ic ic * ic * ie ic J( 1c ic * * * ic ic * Jc * ic -k * i< "k -k * -k 1c * ic -k it * ic -k 1c ic * ic -k ic "k * ic -k 

c * 

C SOLVE FOR THERMAL CONDUCTIVITY AND VISCOSITY OF GAS MIX * 

C * 

C*********************************************************** 

C 

VISC=VISC1/ ( 1 . 0+PHI12* ( 1 . 0-GASMIX) /GASMIX) + 

+ VISC2/ (1. 0+PHI21* (GASMIX/ (1. 0-GASMIX) ) ) 

VISCW=VISCW1/ ( 1 . 0+PHI12* ( 1 . 0-GASMIX) /GASMIX) + 

+ VISCW2/ (1. 0+PHI21* (GASMIX/ (1. 0-GASMIX) ) ) 

KGAS=K1/ (1.0+PSI12* (1. 0-GASMIX) /GASMIX) + 

+ K2/ (1. 0+PSI21* (GASMIX/ (1. 0-GASMIX) ) ) 

C 

c*********************************************************** 

C * 
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SOLVE FOR SPECIFIC HEAT OF GAS MIX 


SOLVE FOR SPECIFIC HEAT OF GAS MIX * 

* 

*********************************************************** 


MW=GASMIX*MW1+ ( 1 . O-GASMIX) *MW2 
CPGAS=GAMMA*RUNIV/ ( (GAMMA-1 .0) *MW) 


C* ********************************************************** 


c * 

C .SOLVE FOR CONVECTION COEFFICIENT * 
C USING THE SIEDER-TATE CORRELATION FOR TURBULENT FLOW * 
C AND THE NUSSELT CORRELATION FOR LAMINAR FLOW * 
C FIND THE NUMBER OF TRANSFER UNITS * 
C * 


Q*********************************************************** 

c 

PR=CPGAS*VISC/ ( KGAS ) 

RE=2 . 0*GASMAS/ ( PI*RICAN*VISC) 

C 

IF (RE .LT. 10000.0) THEN 
F=64 . 0/RE 
NU=3 . 66 

ELSE 

F-(1.82*ALOG10 (RE ) -1 . 64 ) ** ( -2 ) 

NU=0 .027 *RE* *0 . 8 *PR**0 .333333* ( VISC/VISCW) **0.14 
ENDIF 
C 

H=KGAS*NU/ (2 . 0*RICAN) 

NTU=H*2 . 0*PI*RICAN*DZ/ (GASMAS*CPGAS ) 

C 

RETURN 

END 
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C * ************************* * *********************************** * 
C * 

C SUBROUTINE CONDUCT * 

C * 

C PROGRAM DEVELOPS FINITE DIFFERENCE MATRIX * 

C FOR CONDUCTION ANALYSIS IN CYLINDRICAL COORDINATES * 

C * 

c*************************************************************** 

c*************************************************************** 

c * 

C REFERENCE: Pantankar, S. V., "Numerical Heat Transfer and * 

C Fluid Flow", McGraw Hill, 1980, pp. 64-74. * 

C * 

q* ************************************************************* * 

c 

SUBROUTINE CONDUCT (BORL, SORL, QINT) 

C 

REAL K , KE , KLIQ , KMET , KN , KP , KS , KSOL , KW , LENGTH , NTU 
C 

COMMON/ BLOKA/ CVLIQ , CVMET , CVSOL , DENLIQ , DENMET , DENSOL , HFL, 
+KLIQ , KMET , KSOL , LENGTH , PI , RI , RICAN , RO , ROCAN , 

+RMELT(15, 10, 10) , ROMELT ( 15 , 10 , 10 ) , VOID, WALLS 
COMMON/ BLOKC/ DR , DTIME , DPHI , DZ , HEAT , QAZ , QZ , TO ( 15 , 10 , 10 ) 
COMMON/ BLOKD/CPGAS , DELTDZ , DELTP , GASMAS , H , NTU , TGASIN , TGASP , 
, THMELT 

COMMON/ BLOKE/ A (15) , B ( 15 ) , C ( 15 ) , D ( 15 ) , T ( 15 , 10 , 10 ) 

COMMON/ BLOKF/ I , J , L , NI , NJ , NL 
C 

DO 50 1=1, NI 
C 

IF (J .EQ. 1) THEN 
AE=0 . 0 
TE=0 . 0 
ELSE 

TE=T ( I , J-l , L) 

ENDIF 


IF (J .EQ. NJ) THEN 
AW=0 . 0 
TW=0 . 0 
ELSE 

TW=T(I, J+l, L) 
ENDIF 


C 

C 


C 

C SET NODAL CONSTANTS FOR INSIDE RING 
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IF (I .EQ. 1) THEN 
R=RICAN 
RN=RI 

VOID= ( 1 . O-DENLIQ/DENSOL) * (1. O-RMELT ( 2 , J , L) ) 

CALL PROPERTY (CV, DEN , 2 , J , KN , L) 

AA=0 . 5* (RN*RN-R*R) *DPHI*KMET/DZ 
AB=0 . 5* (RN*RN-R*R) *DPHI*KMET/DZ 

AN=KMET*KN*RN*DPHI*DZ/ (0. 5*KMET*DR+KN* (RI-RICAN) ) 

AS=0 . 0 

IF ( J . GT. 1) 

& AE=KMET* (RI-RICAN) *DZ/ (0. 5* (RI+RICAN) *DPHI) 

IF ( J .LT. NJ) 

& AW=KMET* (RI-RICAN) *DZ/ (0. 5* (RI+RICAN) *DPHI) 
AP0=DENMET*CVMET*0 . 5* (RI*RI-RICAN*RICAN) *DPHI*DZ/DTIME 
TCIJ=TGASIN-THMELT 

HI J=GASMAS *CPGAS * ( 1 . 0-EXP (-NTU) ) /FLOAT (2*NJ) 

QIJ=QINT/ (2 . 0 * PI *R* LENGTH) 

C 

0 ************************************************************** 
C * 

C SET NODAL CONSTANTS FOR OUTSIDE RING * 

C * 

C************************************************************** 

c 

ELSEIF (I .EQ. NI) THEN 
R=ROCAN 
RS=RO 

VOID= ( 1 . O-DENLIQ/DENSOL) * ( 1 . O-RMELT ( NI-1 , J , L) ) 

CALL PROPERTY ( CV , DEN , 1-1 , J , KS , L) 

AA=0 . 5* ( R*R— RS *RS ) *DPHI*KMET/DZ 
AB=0 . 5* (R*R-RS*RS) *DPHI*KMET/DZ 

AS=KMET*KS*RS*DPHI*DZ/ (0. 5*KMET*DR+KS* (ROCAN-RO) ) 

AN=0 . 0 

IF ( J .GT. 1) 

& AE=KMET* (ROCAN-RO) *DZ/ ( ( . 5* (ROCAN-RO) +RO) *DPHI) 

IF ( J .LT.NJ) 

& AW=KMET* (ROCAN-RO) *DZ/ ( ( . 5* (ROCAN-RO) +RO) *DPHI) 
AP0=DENMET*CVMET*0 . 5* (ROCAN*ROCAN-RO*RO) *DPHI*DZ/DTIME 
TCIJ=0 . 0 
HIJ=0 . 0 
QIJ=QAZ 


C * 
C SET NODAL CONSTANTS FOR PHASE CHANGE MATERIAL * 
C * 

C************************************************************** 
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C 

ELSE 

R=RI+ ( FLOAT ( I ) -1 - 5 ) *DR 
RN=R+0 . 5* DR 
RS=R-0. 5*DR 

VOID= ( 1 . O-DENLIQ/DENSOL) * ( 1 . O-RMELT ( I , J , L) ) 

CALL PROPERTY (CV, DEN, I , J, K, L) 

AA=0 . 5* (RN*RN-RS*RS ) *DPHI*K/DZ 
AB=0 . 5* (RN*RN-RS*RS) *DPHI*K/DZ 
C 

IF (I .EQ. 2) THEN 

AS=KMET*K*RI*DPHI*DZ/ (0. 5*KMET*DR+K* (RI-RICAN) ) 
TS=T ( 2 , J , L) +DR*AS* (T ( 1 , J , L) -T ( 2 , J , L) ) / 

/ (2 . 0*K*RI*DPHI ) + 

+ 0.125*(T(2,J,L)-2. 0*T ( 3 , J , L) +T ( 4 , J , L) ) 

ELSE 

CALL PROPERTY ( CV, DEN , 1-1 , J , KS , L) 

KP=2 . 0*K*KS/ ( K+KS ) 

AS=KP*RS * DPHI * DZ/ DR 
TS=0.5*(T(I,J,L)+T(I-1,J,L) )+ 

+ 0.125*(T(I+1,J,L)-2.0*T(I,J,L)+T(I-1,J,L) ) 

ENDIF 
C 

IF (I .EQ. NI-1) THEN 

AN=KMET*K*RO*DPHI*DZ/ (0. 5*KMET*DR+K* (ROCAN-RO) ) 
TN=T (NI-1 , J , L) +DR*AN* (T (NI , J , L) -T(NI-1, J,L) )/ 

/ ( 2 . *K*RO*DPHI ) +0 . 125* (T (NI-1 , J , L) - 

2 . 0*T (NI-2 , J , L) +T (NI-3 , J , L) ) 

ELSE 

CALL PROPERTY (CV, DEN, 1+1 , J , KN, L) 

KP=2 . 0*K*KN/ (K+KN) 

AN=KP*RN*DPHI*DZ/DR 
TN=0.5*(T(I,J,L)+T(I+1,J,L) )+ 

+ 0. 125*(T(I+1, J,L) -2.0*T(I,J,L) +T ( 1-1 , J , L) ) 

ENDIF 


C 

C * 
C Calculate the melt front movement * 
C * ********** 


C 

IF (TS*TN .GE. 0.0) THEN 

OVRRLX = BORL*ABS (T ( I , J , L) ) 

DMELT = OVRRLX* 0. 5* (TN+TS) *CV/HFL 
IF (TS .GE. 0.0) THEN 

DMELT=AMIN1 ( DMELT , 1 . O-RMELT ( I , J , L) ) 
ELSE 

DMELT=AMAX1 ( DMELT , 0 . O-RMELT ( I , J , L) ) 
ENDIF 
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ELSE 

IF (T ( I , J , L) *TN .GE. 0.0) THEN 
IF (TN .GE. 0.0) THEN 

DMELT= ( 0 . 5*TN/ (TN-T ( I , J , L) ) -RMELT ( I , J , L) ) *SORL 
ELSE 

DMELT= ((0.5 *TN-T ( I , J , L) )/ 

/ (TN-T ( I , J , L) ) -RMELT ( I , J , L) ) *SORL 

ENDIF 
ELSE 

IF (TN .LE. 0.0) THEN 

DMELT= ( . 5*TS/ (TS-T(I, J,L) ) -RMELT ( I , J , L) ) *SORL 
ELSE 

DMELT= ((0.5 *TS-T ( I , J , L) ) / 

/ (TS-T ( I , J , L) ) -RMELT ( I , J , L) ) *SORL 

ENDIF 
ENDIF 
ENDIF 
C 

C************************************************************** 

c * 

C Check melt front movement to ensure that it meets * 

C thermodynamic limitations. * 

C * 

c************************************************************** 

c 

IF (DMELT .GE. 0.0) THEN 
DO 40 11=1+1 , NI-1 

IF (RMELT (II, J, L) .LT. 0.99) DMELT=0 . 0 
40 CONTINUE 

ENDIF 
C 

IF (I .EQ. 2) THEN 

IF (DMELT .LT. 0.0) THEN 

IF( (RMELT (I+1,J, L) .LE. .01) .OR. 

& (T ( 1-1 , J , L) .LE. 0.0)) 

& RMELT ( I , J , L) =RMELT ( I , J , L) +DMELT 

ELSE 

IF ( (RMELT (I+1,J, L) . GE . 0.99) .OR. 

& (T ( 1-1 , J , L) .GE. 0.0)) 

& RMELT ( I , J , L) =RMELT ( I , J , L) +DMELT 

ENDIF 

ELSEIF ( I .EQ. NI-1 ) THEN 
IF (DMELT .GE. 0.0) THEN 

IF ( (T (NI , J, L) .GE. 0.0) .OR. 

& (QAZ .GE. 0.0)) 

& RMELT ( I , J , L) =RMELT ( I , J , L) +DMELT 

ELSE 

IF ( (T (NI , J , L) .LE. 0.0) .OR. 

& . (RMELT (NI-2,J,L) .LE. .01)) 

& RMELT ( I , J , L) =RMELT ( I , J , L) +DMELT 



U? ft? ft*> 


ENDIF 
ELSE 

IF ( DMELT .LT. 0.0) THEN 

IF( (RMELT(I-1, J, L) .LE. .01) .OR. 

(RMELT(I+1, J, L) .LE. .01)) 

RMELT ( I , J , L) =RMELT ( I , J , L) + DMELT 

ELSE 

IF ( (RMELT (I-1,J,L) . GE . .99) .OR. 

(RMELT ( 1+1 , J , L) .GE. .99)) 

RMELT ( I , J , L) =RMELT ( I , J , L) + DMELT 

ENDIF 
ENDIF 

RMELT ( I , J , L) =AMAX1 (0.0, RMELT ( I , J , L) ) 

RMELT ( I , J , L) =AMIN1 (1.0, RMELT ( I , J , L) ) 

IF (J .GT. 1) THEN 

CALL PROPERTY ( CV, DEN , I , J-l , KE , L) 

KP=2 . 0*K*KE/ (K+KE) 

AE=KP*DR*DZ/ (R*DPHI ) 

ENDIF 

IF ( J .LT. NJ) THEN 

CALL PROPERTY ( CV, DEN , I , J+l , KW, L) 

KP=2 . 0*K*KW/ (K+KW) 

AW=KP*DR*DZ/ (R*DPHI ) 

ENDIF 

APO=DEN*CV*R*DPHI*DR*DZ/DTIME 
TCI J=0 . 0 
HI J =0 . 0 

PCMMASS=R*DPHI*DR*DENLIQ* (1. 0-WALLS) *DZ 
QI J= (ROMELT ( I , J , L) -RMELT ( I , J , L) ) * PCMMASS *HFL/ DTI ME 
ENDIF 

IF (L .EQ. 1) THEN 
AA=0 . 0 
TA=0 . 0 
ELSE 

TA=T ( I , J , L-l) 

ENDIF 

IF (L .EQ. NL) THEN 
AB=0 . 0 
TB=0 . 0 
ELSE 

TB=T ( I , J , L+l) 

ENDIF 


B ( I ) =AA+AB+AN+AS+AE+AW+APO+HI J 
C ( I ) =-AN 

D ( I ) =APO*TO ( I , J , L) +HIJ *TCIJ +QIJ+AE*TE+AW*TW+AA*TA+AB*TB 

CONTINUE 

RETURN 

END 
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C************************************************************** 
C*****************************************************' ********* 

c * 

C SUBROUTINE PROPERTY * 

C * 

C This program calculated the effective thermal * 

C properties of the solid & liquid phases of the storage media* 
C * 

c************************************************************** 
c************************************************************** 
c*** ******************************************** *************** 

c * 

C REFERENCES Y. S. Touloukian et.al. "Thermophysical * 

C Properties of Matter", Vol. 2, 1970. * 

C * 

C R. P. Wichner, ORNL "Boundary Conditions and Thermophysical * 
C Properties of LiF-CaF2 " , Presented to NASA Lewis Research * 
C Center, June 10, 1987. * 

C • * 

C R. I. Vachon, et.al, "Thermal Conductivity of Heterogeneous * 
C Mixtures & Lunar Soils", NAS8-26579, October 1973. * 

C * 

Q*** **************************************************** ******* 

C 

SUBROUTINE PROPERTY ( CV , DEN , I , J , K, L) 

C 

COMMON/ BLOKA/CVLIQ, CVMET, CVSOL, DENLIQ, DENMET, DENSOL, HFL, 
+KLIQ , KMET , KSOL , LENGTH , PI , RI , RICAN , RO , ROC AN , 

,RMELT(15, 10, 10) , ROMELT ( 15 , 10 , 10 ) , VOID, WALLS 
C 

REAL K , KLIQ , KMET , KP , KPSOL , KPURE , KSOL , LENGTH 


C******* ************************************************* ****** 

c * 

C FIND THE THERMAL PROPERTIES OF THE PURE MATERIALS * 

C * 

C CALCULATE EFFECT OF VOID ON PROPERTIES OF SOLID * 

C USE CORRELATION BY MEREDITH & TOBIAS TO CORRECT CONDUCTIVITY* 
C * 

c************************************************************** 

c 

KPSOL=KSOL* (2 . 0* (1. 0-VOID) +0. 6 13 5* VOID* *2 .333-1. 6*VOID** 

& 3 . 333)/ (2 . 0+VOID+0. 6135*VOID**2 . 333-0 . 6975*VOID**3 .333) 
DENPSOL= DENLIQ 
C 

KPURE=1 . 0/ (RMELT ( I , J , L) /KLIQ+ ( 1 . 0-RMELT ( I , J , L) ) /KPSOL) 
CVPURE=RMELT ( I , J , L) *CVLIQ+ ( 1 . 0-RMELT ( I , J , L) ) *CVSOL 
DENPURE=RMELT ( I , J , L) *DENLIQ+ ( 1 . 0-RMELT ( I , J , L) ) *DENPS0L 


o o o o o o o 
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************************************************************** 

* 

PUT PURE STORAGE MATERIAL IN PARALLEL WITH CAN MATERIAL * 

* 

************************************************************** 

K=KPURE* ( 1 . O-WALLS) +KMET*WALLS 
DEN=DENPURE* ( 1 . 0 -WALLS ) +WALLS*DENMET 

CV= (CVPURE*DENPURE* ( 1 . 0 -WALLS ) +WALLS*CVMET*DENMET) /DEN 
C 

RETURN 

END 



ooooo oooooo 


c************************************************************** 


c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE SOLVER 

This program solves a tridiagonal matrix using 
an elimination method 

[ A] T ( 1-1) + [ B] T ( I ) + [C]T(I+1) = [D] 


* 

* 

* 

* 

* 

* 

* 

* 


c************************************************************** 

c 

SUBROUTINE SOLVER 

COMMON/BLOKE/A ( 15 ) , B ( 15 ) , C ( 15) , D ( 15 ) , T ( 15 , 10 , 10 ) 

COMMON/ BLOKF/ I , J , L , NI , NJ , NL 


PROCEED WITH GAUSSIAN ELIMINATION TO DEVELOP A * 

UPPER DIAGONAL MATRIX * 

* 


DO 1, 1=1 , NI— 1 

R=A(I+1)/B(I) 

B ( 1+1 ) =B ( 1+1 ) -R*C ( I ) 
D ( 1+1 ) =D ( 1+1 ) -R*D ( I ) 
1 CONTINUE 


BEGIN BACK SUBSTITUTION * 

* 


C 

T (NI , J , L) =D (NI ) / B (NI ) 

C 

DO 2, I=NI-1 ,1,-1 

T(I,J,L)=(D(I)-C(I)*T(I+1,J,L) )/B(I) 
2 CONTINUE 
RETURN 
END 


APPENDIX B 


PROGRAM FOR THE THERMAL ANALYSIS 
OF THE 

CAVITY RECEIVER 
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C PROGRAM SOLDTESS * 

C * 

C SOLAR DYNAMIC TRANSIENT ENERGY STORAGE SIMULATION * 

C * 

C This program is the main routine for the TESS code. * 

C Its purpose is to read & print data and to call * 

C calculational subroutines so as to complete the thermal * 

C analysis. Checks for convergence are made here. * 


c******************************************************************* 

c 

C NOMENCLATURE 

C 

•C System dimensions 

C CS = Coefficient for Simpson 1/3 Integration 

C DR = Diff. element dimension in radial direction, m 

C DPHI = Diff. element dim. in azimuthal direction, rad 

C DMAX = Maximum move in melt position in one iteration 

C DMELT = Fractional movement of melt position in element 

C DTIME = Time increment, sec 

C DZ = Diff. element dimension in axial directrion, m 

C GASMAS = Coolant flow rate, kg/sec 

C LENGTH = Axial length of stack, m 

C OVRRLX = Over relaxiation factor on melt movement 

C POS = Axial position of current node, m 

C R = Radius of element (I, J) , m 

C RI = Outside radius of inner can ring, m 

C RICAN = Inside radius of inner can ring, m 

C RMELT = Fraction of PCM in element melted 

C RN = Radius at outside edge of element (I, J) , m 

C RO = Inside radius of outer can ring, m 

C ROCAN = Outside radius of outer can ring, m 

C RS = Radius at inside edge of node (I, J) , m 

C ROMELT = Previous fraction of PCM melted 

C T = Projected relative temperature of node 

C TCIJ = Absolute coolant temperature of gas at pos., K 

C TGAS = Relative coolant temperature at axial position 

C TGASIN = Absolute coolant temp, entering axial position, K 

C TGASP = Absolute coolant temperature at mid DZ , K 

C THGAS = Absolute coolant temperature at inlet, K 

C THETA = Absolute temperature of node (I, J) , K 

C THETAG = Absolute temperature of gas, K 

C THMELT = Absolute melting temperature of storage media, K 

C TIME = Time into orbit run, min 

C TIMER = Time of orbit, min 
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TO = Previous relative temperature of node (I,J) 

TP = Projected relative temperature of node (I,J) 

WALLS = Ratio of wall conduction area to storage area 

Material properties 

CPGAS = Coolant specific heat, kJ/kg K 
CVLIQ = Liquid specific heat, kJ/kg K 
CVMET = Ring specific heat, kJ/kg K 

CVPURE = Specific heat of PCM (solid & liquid) kJ/kg K 
CVSOL = Solid specific heat, kJ/kg K 
DENLIQ = Liquid density, kg/m**3 
DENMET = Ring density, kg/m**3 

DENPURE = Density of PCM (solid & liquid) kg/m**3 
DENSOL = Solid density, kg/m**3 
F = Smooth tube adiabatic friction factor 
GASMIX = Fraction of gas 1 in coolant mix 
H = Convection coefficient of gas W/m**2 K 
HFL = Heat of fusion, kJ/kg 

K = Eff. thermal conductivity of node (I,J), W/m K 
KE = Thermal conductivity of node (I,J+1) , W/m K 

KN = Thermal conductivity of node (I+1,J) , W/m K 

KP = Thermal conductivity solid with void, W/m K 
KS = Thermal conductivity of node (I-1,J) , W/m K 

KW = Thermal conductivity of node (I,J-1), W/m K 

K_1 = Thermal conductivity constant (See GASPROP) 

K_2 = Thermal conductivity constant (See GASPROP) 

K_3 = Thermal conductivity constant (See GASPROP) 

K_4 = Thermal conductivity constant (See GASPROP) 

KLIQ = Liquid thermal conductivity, W/m K 
KMET = Ring thermal conductivity, W/m K 

KPURE = Eff. thermal cond. of PCM (solid & liq. ) W/m K 
KSOL = Solid thermal conductivity, W/m K 
MW_ = Molecular weight of gas component _ 

VISC_1 = Viscosity constant (see GASPROP) 

VISC_2 = Viscosity constant (see GASPROP) 

VISC_3 = Viscosity constant (see GASPROP) 

VISC_4 = Viscosity constant (see GASPROP) 

VOID = Void fraction in element 

Indicies 

I = Counter in radial direction 
ICOUNT = Loop counter 

I FLAG = Flag to indicate insufficient convergence 

IMAX = Maximum number of loops 

J = Counter in azimuthal direction 

L = Counter in axial direction 

LIO = Counter on output file 

NORBIT = Number of orbit cycles 

NI = Number of radial nodes 

NJ = Number of azimuthal nodes 


B3 


C NL = Number of axial nodes 

C NTIME = Number of time increments per orbit 

C 

C Matrix coefficients 

C A = Conduction matrix coefficient on 1-1 node 

C B = Conduction matrix coefficient on I node 

C C = Conduction matrix coefficient on 1+1 node 

C D = Conduction matrix coefficient on constant vector 

C 

C Thermal resistances & capacitances 

C AE = Thermal conductance of J-l node, W/K 

C AN = Thermal conductance of I+i node, W/K 

C AS = Thermal conductance of 1-1 node, W/K 

C AW = Thermal conductance of J+l node, W/K 

C APO = Thermal capacitance of node (I,J) , kJ/m K 

C 

C Heat fluxes 

C QAZ = Outside surface heat flux, W/m**2 

C QDIR = Direct solar heat flux, W/m**2 

C QGASP = Total heat added to gas at position DZ/2, kJ 

C QGASDZ = Total heat added to gas at post ion DZ , kJ 

C QIJ = Radiation to element (I,J) , W/m* ’ : 2 

C QIDIR = Indirect solar heat flux, W/m**2 

C QINTR = Heat input on internal surface, W/m**2 

C QZ = Integrated avg. axial multiplier on heat flux 

C QZZ = Local axial multiplier on heat flux 

C QZ_ = Constant on fit of QZ 

C 

Q* ****************************************************************** 

C * 

PROGRAM SOLDTESS 


C 

REAL KLIQ , KMET , KSOL, LENGTH , MIN ( 11 ) , NTU 
DIMENSION THETA ( 10 ), TP (15, 10) 

C 

COMMON/ BLOKA/CVLIQ, CVMET , CVSOL, DENLIQ , DENMET , DENSOL, HFL, 

, KLIQ, KMET, KSOL, LENGTH, RMELT ( 25, 20) , ROMELT ( 25 , 2 0 ) ,VOID 
COMMON/ BLOKC / DT IME , DX , DY , TO ( 2 5 , 2 0 ) 

COMMON/ BLOKE/ A (25) ,B(25) ,C(25) ,D(25) ,Q(10) ,QNA(20) ,T(25,20) 
COMMON/ BLOKF/ I , J , N I , N J 
DATA U/400.0/ 

C 

£****************************************************************** 
C * 

c OPEN INPUT AND OUTPUT FILES * 

C * 

Q****************************************************************** 

C 

OPEN ( 8 , FILE= ' INPUT . TXT ' ) 

OPEN ( 7 , FILE= ' OUTPUT . TXT ' ) 
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11=1 

READ (8,2100, ERR=93 ) IMAX, NORBIT , NI , NJ , NTIME 
READ (8,2200, ERR=9 4 ) H , HI , W , BORL , SORL 

READ (8,2200, ERR=9 5 ) CVLIQ , DENLIQ , KLIQ , LENGTH , THMELT 
READ (8,2200, ERR=9 6 ) CVSOL , DENSOL , KSOL , ERRMAX , HFL 
READ (8,2200, ERR=9 7 ) CVMET , DENMET , KMET , GASMIX , TIMER 
1 READ (8,2200, ERR=99 , END=5 ) MIN ( II ) ,Q(II) , DUM2 , DUM3 , DUM4 
11=11+1 
MIN (II) =1 . E9 
GO TO 1 
C 

Q********************************************* -k ******************** 


c * 

C INITIALIZE LOOP COUNTERS AND ERROR FLAGS * 

C * 

c****************************************************************** 
C 

5 ITEMP=0 

DX=H1/FL0AT ( NJ ) 

NI=INT (H/DX+0 . 5 ) 

DXX=H-DX*FLOAT (NI-1 ) 

DY=W/FLOAT(2*NJ) 

DTIME=6 0 . 0 *TIMER/ FLOAT (NTIME ) 

C 


DO 10 11=1, NI 
DO 10 JJ=1 , NJ 
T (II , JJ) =-l . 0 
TP (II , JJ) =T ( II , JJ) 
RMELT ( II , JJ) =0 . 0 
10 CONTINUE 

T1AVG=T (1,1) 




C * 

C BEGIN CALCULATIONAL LOOPS ON LENGTH, TIME, ANGLE * 

C * 

C 

DO 100 NO=l, NORBIT 

NPUT=0 

TIME=0 . 0 

DO 100 N=l, NTIME 

C 

IF (MIN (NPUT+1) .LE. TIME+0 . IE-6) NPUT=NPUT+1 

C 


DO 20 JJ=1 , NJ 

QNA(JJ) =Q (NPUT) /FLOAT (NJ) 
DO 20 11=1, NI 

ROMELT ( II , JJ) = RMELT ( I I , J J ) 
TO ( II , JJ) =T ( II , JJ) 



20 CONTINUE 

C 

ICOUNT=0 
30 IFLAG=0 

ICOUNT=ICOUNT+l 

C 

CALL SODIUM ( NPUT , THMELT , TKNA , TNA , U) 

C 

DO 40 J=1,NJ 

CALL CONDUCT (BORL, DXX , SORL, TNA , U) 

CALL SOLVER 
40 CONTINUE 

C 

DO 50 11=1 , NI 
DO 50 JJ=1 , NJ 

ERROR=ABS ( T ( I I , J J ) -TP (II, JJ) ) /THMELT 
IF (ERROR .GT. ERRMAX) IFLAG=2 
IF (ERROR .GT. ERRMAX) ITEMP=ITEMP+1 
TP ( II , JJ) =T (II , JJ) 

50 CONTINUE 

C 

I F ( I COUNT . GT . IMAX ) THEN 
WRITE (7 , 1000) ITEMP 
IFLAG=0 
ENDIF 

C 

IF ( I FLAG .GE. l)GO TO 30 
TIME=TIME+DTIME/60 . 0 
IF (NO .LT. NORBIT) GO TO 100 
C 

WRITE(7, 1050) 

WRITE (7, 1250) 

WRITE (7 , 1100) TIME, TKNA 
DO 70 11=1, NI 

DO 60 JJ=1 , NJ 

THETA ( J J ) =T ( I I - ( J J- 1 ) , J J ) +THMELT 
60 CONTINUE 

JJJ=AMIN0 ( II , NJ ) 

WRITE(7, 1200) (THETA (JJ) , JJ=1, JJJ) 

70 CONTINUE 

C 

WRITE(7, 1050) 

WRITE ( 7 , 1250) 

DO 80 1=1, NI 

JJJ=AMINO(NJ, I) 

WRITE(7, 1350) (RMELT ( I- ( JJ-1) , JJ) , JJ=1, JJJ) 
80 CONTINUE 

C 

100 CONTINUE 
C 



o o o 
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c*************************************************************** 

* 

INPUT AND OUTPUT FORMAT STATEMENTS * 

* 

C** ************************************************************* 

c * 

1000 FORMAT (IX, 'FAILED TO CONVERGE AFTER ',16,' ITERATIONS.') 

1050 FORMAT (IX, ' ') 

1100 FORMAT ( IX, ' Time = ', 1F6.1,' min. Sodium Temp. = ', 

, 1F6 . 1 , ' K.') 

1200 FORMAT ( 10F7 . 1) 

+ ' ==================== ' ) 

1300 FORMAT (IX, 'Avg. Gas Temp. = ',F6.1,'K Out. Gas temp . ' , F6 . 1 , ' K H = 
+ ' ,F7.4, 'KW/m2-K') 

1350 FORMAT (10F7 .4) 

1400 FORMAT (IX, 'INPUT FORMAT ERROR IN 1ST LINE') 

1401 FORMAT (IX, 'INPUT FORMAT ERROR IN 2ND LINE') 

1402 FORMAT (IX, 'INPUT FORMAT ERROR IN 3RD LINE') 

1403 FORMAT (IX, 'INPUT FORMAT ERROR IN 4TH LINE') 

1404 FORMAT (IX, 'INPUT FORMAT ERROR IN 5TH LINE') 

1405 FORMAT (IX, 'INPUT FORMAT ERROR IN 6TH LINE') 

1406 FORMAT (IX, 'INPUT FORMAT ERROR IN TRANSIENT SECTION') 

2100 FORMAT (16, 4110) 

2200 FORMAT (5F1 0.4) 

STOP 

93 WRITE(*, 1400) 

STOP 

94 WRITE(*, 1401) 

STOP 

95 WRITE(*, 1402) 

STOP 

96 WRITE(*, 1403) 

STOP 

97 WRITE (*, 1404) 

STOP 

98 WRITE(*, 1405) 

STOP 

99 WRITE(*, 1406) 

STOP 

END 


APPENDIX C 


PROGRAM FOR THE THERMAL ANALYSIS 
OF THE 

PEBBLE BED RECEIVER 



Cl 


c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


c 


c 


PROGRAM MAIN 

C********************************************************* 


C * 

C TEST COEFF * 

C * 

C * 


C ******************************************************* ** 
C* ******************************************************** 

c 

DIMENSION TL(15, 10, 20) ,RL(15, 10, 20) 

COMMON/ BLOKF/ I , J , L , NI , NJ , NL , ITIME , MODE 

COMMON/BLOKA/RMELT (15,10), KLIQ , KSOL , KG AS , KMET , VOI D , EPS , 

, DRM , RO , DENLIQ , CVLIQ , CVSOL , DENMET , CVMET , HFL, 

, S APUV , PI , PR , REP , H , CVPURE , DENPURE , DENSOL, 

, ROMELT ( 15 , 10) , RPMELT ( 15 , 10) 

COMMON/ BLOKB/DZ , DTI ME , DR, QAZ , QZ , HEAT 

COMMON/ BLOKD/TGAS (50) ,TO(15, 10) ,T(15, 10) , AREA , GASMAS , 

, TGASO ( 50 ) , G ASMIX , TG AS IN , DP , VF , DPHI , CPGAS , 

, DENGAS , TP ( 15 , 10 ) ,TGASP(50) , TGAST ( 50 , 100) 

COMMON/ BLO KG/ KE ( 15 , 10 ) ,KW(15, 10) , KN , KS , HWGAS ( 10) ,HWSOL 
COMMON/BLOKJ/A( 15) , B ( 15 ) , C ( 15) , D ( 15 ) 

COMMON/ BLO KO/KPSOL, KPURE , K ( 3 , 3 ) 

COMMON/ BLOKX/ FCTR 

COMMON/BLOKS/Q1DIR, Q1IN, Q2DIR, Q2IN 

REAL KPSOL , KSOL , KPURE , DFT , DF , KSTAG , KCOND , KRAD , KRO , FS , KWKF 

REAL KFSTAG , KFF , KWALL , KWRAD , K , KLIQ , KGAS , KMET , LENGTH 

REAL KN , KS , KE , KW 

INTEGER DTIME 

OPEN ( 8 , FI LE= ' INPDATA . TXT ' ) 

READ ( 8 , * ) NI , NJ , NL, NTIME , DTIME , NORBIT 

READ ( 8 , * ) DP , DRM , RO , LENGTH , VF , S APUV , S APUVW 

READ ( 8 , * ) KLIQ , CVLIQ , DENLIQ , THMELT 

READ ( 8, *) KSOL, CVSOL, DENSOL, HFL 

READ ( 8 , * ) KMET , CVMET , DENMET , EPS 

READ ( 8 ,*) TGAS IN , GASMAS , GASMIX 

READ (8, *) FCTR, ERRMAX, Q1DIR, Q1IN, Q2DIR, Q2IN, IMAX 

PI=3 . 14159 
DPHI=3 . 14159/NJ 
AREA=PI*RO**2 
DR=RO/ (NI-1) 

DZ=LENGTH/NL 

OPEN ( 5 , FILE= 'TEMP . TXT ' ) 

DO 12 IORBIT=l, NORBIT 
DO 12 MODE=l , 2 

IF (MODE .EQ. 2) NTIME=2160 
IF (MODE .EQ. 1) NTIME=3 2 40 



DO 12 L=1,NL 

DO 12 ITIME=1 , NTIME/DTIME 

IF ( ( ITIME .EQ. 1) .AND. (MODE .EQ. 1) 
& .AND. (IORBIT .EQ. 1)) THEN 

TGAS(L)=TGASIN 
TGASO ( L) =TGASIN 
DO 6 11=1, NI 
DO 6 JJ=1,NJ 

T ( II , JJ) =1218.0 
TO (II , JJ) =1218.0 
TP ( II , JJ) =1218.0 
ROMELT (II , JJ) =0 . 0 
RMELT ( II , JJ) =0 . 0 

6 CONTINUE 

ELSEIF (ITIME .EQ. 1) THEN 

TGASO ( L) =TGAST ( L, NTIME/DTIME) 
TGAS(L)=TGASO(L) 

TGASP(L) =TGASO(L) 

DO 7 11=1, NI 
DO 7 JJ=1,NJ 

TO ( II , JJ ) =TL ( II , JJ , L) 

T ( II , JJ) =TO (II , JJ) 

TP ( II , JJ) =TO ( II , JJ) 

ROMELT ( II , JJ) =RL( II , JJ , L) 

RMELT ( I I , J J ) = ROMELT ( I I , J J ) 

7 CONTINUE 


400 

ENDIF 
JKJ=0 
HEAT=0 . 0 

500 

DO 11 J=1 , NJ 

CALL TGAS2 (KGAS , H, PR, REP) 
DO 10 1=1, NI 

10 

CALL TCON 
CALL TCOEFF 
CONTINUE 

11 

CALL TSOLVER 
CALL TMELT ( THMELT ) 
CONTINUE 

2200 

CALL TGASTEMP 

FORMAT (IX, ' TGAS ( ' ,11, ') = ' , F7 


DO 30 J=1,NJ 


DO 30 1=1, NI 

IF ( ABS ( T ( I , J ) -TP ( I , J ) ) .GE. ERRMAX) THEN 
DO 40 JJ=1,NJ 
DO 40 11=1, NI 

TP ( II , JJ) =T ( II , JJ) 

40 CONTINUE 

JKJ=JKJ+1 

WRITE (*, 2216) IORBIT, MODE, L, ITIME, JKJ 
FORMAT ('+',413 , 5X, 14) 


2216 
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IF ( JKJ .GT. IMAX) THEN 
WRITE (5,7000) JKJ 

7000 FORMAT (IX, 'FAILED TO CONVERGE AFTER', 14,' 

& ITERATIONS') 

GO TO 700 
ENDIF 
GO TO 400 
ENDIF 

30 CONTINUE 

SUMT2=0 . 0 
SUMR2=0 . 0 
DO 28 1=1, NI 
SUMT1=0. 0 
SUMR1=0 . 0 
DO 29 J= 1 , NJ 

CVPURE=RMELT ( I , J) *CVLIQ+ 

+ ( 1 . 0-RMELT ( I , J) ) *CVSOL 

SUMT1=SUMT1+ (T (I , J) -TO ( I , J) ) *CVPURE 
IF (I .EQ. 1) GO TO 29 
SUMR1=SUMR1+ ( RMELT ( I , J) -ROMELT ( I , J) ) 

29 CONTINUE 

RS= (NI-I ) *RO/ (NI-1 ) 

RN=RS+RO/ (NI-1) 

VOLUME= ( DPHI/2 ) * (RN-RS) * (RN+RS) *DZ 
PRO DT=VO LUME * SUMT 1 
PRODR=VOLUME * SUMR1 
SUMT2=SUMT2+PRODT 
SUMR2=SUMR2+PRODR 
28 CONTINUE 

QLAT=2* ( 1-VF) *DENPURE*SUMR2*HFL/DTIME 
QSENS=2 * ( 1-VF) *DENPURE*SUMT2/DTIME 
IF (L .EQ. 1) THEN 
TG AS B=TG AS I N 
ELSE 

TGASB=TGAST ( L-l , ITIME) 

ENDIF 

QG AS =G ASMAS * C PG AS * (TGAS (L) -TGASB) 


700 


3500 

& 


2000 


TGASO(L) =TGAS (L) 

TGAST ( L , ITIME ) =TGAS ( L) 

WRITE (5,*) 

WRITE (5,*) 

WRITE (5, *) 

WRITE ( 5 , 3500) IORBIT , MODE , ITIME , L 

F0RMAT(1X, 'IIIIII ORBIT # ',11,' IIIIII MODE ',11,' 
IIIIII TIME= ',11, ' IIIIII L= ',11,' IIIIII') 

WRITE (5, *) 

DO 15 1=1, NI 

WRITE(5, 2000) (T(I,J) ,J=NJ,1,-1) 

FORMAT (IX, 10 (F6. 1, IX) ) 
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15 CONTINUE 
WRITE (5, *) 

DO 18 1=1, NI 

WRITE (5, 2100) (RMELT ( I , J) ,J=NJ, 1,-1) 

2100 FORMAT ( IX ,10(F6.3,1X) ) 

18 CONTINUE 

WRITE (5, *) 

WRITE (5,8000) KGAS , H , PR, REP 

8000 FORMAT (IX, ' KGAS= ' , F9 . 7 , 3X , 'H=' , F6 . 4 , 3X, 'PR=',F6.4, 

& 3X , ' REP= ' , F8 . 3 ) 

WRITE ( 5, 2200 )L,TGAS(L) 

HEAT=HEAT*DZ 
WRITE (5, 8001) QZ, HEAT 

8001 FORMAT (IX, 'HEAT INPUT FRAC. = ' , F8 . 5 , 3X , ' HEAT= ' , F8 . 5 ) 
WRITE (5, *) 

WRITE ( 5 ,8002 ) QGAS , QLAT , QSENS 

8002 FORMAT (IX, ' QGAS= ' , F7 . 5 , 5X , ' QLAT= ' , F7 . 5 , 5X , 

& ' QSENS= ' , F7 . 5 ) 

DO 16 1=1, NI 

DO 16 J=1 , NJ 

TO ( I , J ) =T ( I , J ) 

TP (I , J) =T (I , J) 

TL ( I , J , L) =T ( I , J ) 

ROMELT ( I , J) = RMELT ( I , J) 

RL ( I , J , L) =RMELT ( I , J) 

16 CONTINUE 
12 CONTINUE 

5000 END 
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C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c* * ******* ******* * * ****** ********************************** 

q* ******************* ******* ******* ******* ***************** 

c * 

C SUBROUTINE GASPROP * 

c * 

C This program calculates the effective thermal properties* 
C and convective coefficients for the gas coolant. * 

C The gas coolant is assumed to be a mixture of He and Xe.* 
C Properties of these two constituents are included as * 
C curve fits with the corresponding coefficients included * 
C as data in the subroutine. * 

C * 

q* ******************************************************** ■* 
q** ******************************************** ************ 


c 

q** ******************************************* ************* 

c * 

C REFERENCES * 

C M. R. Vanco, "Analytical Comparison of Relative Heat * 
C Transfer Coefficients & Pressure Drops of Inert Gases * 
C and Their Binary Mixtures", Lewis Research Center * 

C * 

C A. J. Chapman "Fundamentals of Heat Transfer" * 

C Macmillan, 1978, p 329. * 

C * 

Q* ******************************************************** * 

c 

c********************************************************** 

c * 

C EQUATIONS * 

C K1 = Kll + K12 * (TGAS - K13)**K14 * 

C K2 = K21 + K2 2 * (TGAS - K23)**K24 * 

C * 

C VISC1 = VISC11 + VISC12* (TGAS - VISC13 ) **VISC14 * 

C VISC2 = VISC21 + VISC22* (TGAS - VISC23 ) **VISC24 * 

C * 

C********************************************************** 

c 

SUBROUTINE TGAS2 ( KGAS , H , PR , REP) 

REAL K, K1 , K2 , Kll , K12 , K13 , K14 , K2 1 , K22 , K23 , K24 , MW , MW1 , MW2 , 

, KGAS 

COMMON/ BLOKD/ TGAS ( 50 ) ,TO(15, 10) ,T(15, 10) , AREA , GASMAS , 

, TGASO ( 50 ) , GASMIX, TGASIN, DP, VF, DPHI , CPGAS , 

, DENGAS , TP ( 15 , 10 ) ,TGASP(50) , TGAST ( 50 , 100) 

COMMON/ BLOKF/ I , J , L, NI , NJ , NL, ITIME , MODE 

COMMON/BLOKG/KE ( 15 , 10) ,KW(15,10) , KN , KS , HWGAS ( 10) , HWSOL 
DATA GAMMA, PI, RUNIV, MW1 , MW2/ 

/l. 66, 3. 14 159, 8. 3 14 34, 4. 0,13 1.0/ 

DATA Kll , K12 , K13 , K14/0 . 1580E-3 , 0 . 9086E-6 ,333.333,0. 794/ 
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C 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 


DATA K2 1 , K22 , K23 , K24/ 0 . 640E-5 , 0 . 0454E-6 ,333.333,0.800/ 

DATA VISC11 , VISC12 , VISC13 , VISC14/ 

/ . 2 158E-4 , . 1022E-6, 333 . 333 , .825/ 

DATA VISC21, VISC22 , VISC23 , VISC24/ 

/0 . 253E-4 , 0 . 2452E-6 ,333.333,0.77/ 

F1(X,Y)=(1. 0+2 .41* (X-Y) * (X-0 . 142*Y ) / (X+Y) **2) 

F2 (X1,X2,X3,X4,X5,X6)= 

= ( ( 1 . 0+ (X1/X2 )**0.5*(X3/X4) **0.25) **2) / 

/(2. 0*2. 0**0. 5*(1. 0+X5/X6) **0.5) 

C 

C* ********************************************************* 

C * 

C FIND THE THERMAL CONDUCTIVITY & VISCOSITY OF PURE GASES * 

C * 

C********************************************************** 

C 

K1=K11+K12* (TGAS (L) -K13) **K14 
K2=K21+K22* (TGAS (L) -K23) **K24 

VISC1=VISC11+VISC12* (TGAS (L) -VISC13 ) **VISC14 
VISC2=VISC2 1+VISC2 2* (TGAS (L) -VISC23) **VISC24 
C 

c ********************************************************** 
C * 

C FIND THE WEIGHTING PARAMETERS FROM REFERENCE * 

C * 

c * ********************************************************* 


C 

PHI12=F2 (VI SCI , VISC2 , MW2 , MW1 , MW1 , MW2 ) 

PHI21=PHI12* (VISC2*MW1) / (VISC1*MW2 ) 

PSI12=F2 ( K1 , K2 , MW1 , MW2 , MW1 , MW2 ) *F1 (MW1,MW2) 

PSI2 1=F2 (K2 , K1 , MW2 , MW1 , MW2 , MW1 ) *F1 (MW2,MW1) 

C 

C* ********************************************************* 

C * 

c SOLVE FOR THE THERMAL CONDUCTIVITY * 

C AND VISCOSITY OF GAS MIX * 

C * 

C* ********************************************************* 


C 

VISC=VISC1/ (1. 0+PHI12 * ( 1 . 0-GASMIX) / GASMIX) + 

/ VISC2/ (1. 0+PHI21* (GASMIX/ (1. 0-GASMIX) ) ) 

KGAS=K1/ ( 1 . 0+PSI12* ( 1 . 0-GASMIX) /GASMIX) + 

/ K2/(1.0+PSI21* (GASMIX/ (1. 0-GASMIX) ) ) 

C* ************************************ ********************* 

c * 

c SOLVE FOR SPECIFIC HEAT OF GASMIX * 

C * 

C* ********************************************************* 


C 

MW=GASMIX*MW1+ ( 1 . 0-GASMIX) *MW2 
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C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


CPGAS=GAMMA*RUNIV/ ( (GAMMA-1.0) *MW) 

C 

C*** ************************* ****************************** 


c * 

C CORRELATION RECOMMENDED FOR GAS FLOW IN A BED OF SPHERES* 
C _ _0 . 575 * 
C (ej) = (ej) = 2.06 (Re) * 
C h m * 
C * 
C j are the Colburn j factors, and e is the void fraction * 
C * 
C REFERENCE * 
C Incropera, Frank P. , "Fundamentals of Heat and Mass * 
C Transfer", 2nd Edition, John Wiley & Sons Inc., 1981. * 
C p 292,352. * 
C * 
C* ******************************************************** * 


c 

PR=CPGAS*VISC/ ( KGAS ) 

REP=GASMAS*DP/ (VISC*AREA) 

H= ( 2 . 06/VF) * ( KGAS/DP) *PR** (1.0/3) *REP**0.425 
DENGAS=11. 22 
C 

RETURN 

END 


i 
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C C 

c c* *************************************************** ******* 

c c*************************************** ******* ************* 


c c * 
c C SUBROUTINE COEFF * 
C C * 
C C This program develops finite difference matrix for * 
C C conduction analysis in cylindrical coordinates. * 
C C * 


c c*********************************************************** 
c c ************************************** . ******************** 
c c 

SUBROUTINE TCOEFF 

COMMON/ BLOKA/RMELT (15,10), KLIQ , KSOL , KGAS , KMET , VOID , EPS , 

, DRM , RO , DENLIQ , CVLIQ , CVSOL , DENMET , CVMET , HFL, 

, SAPUV , PI , PR , REP , H , CVPURE , DENPURE , DENSOL , 

, ROMELT (15,10) , RPMELT (15,10) 

COMMON/ BLOKB/DZ , DTIME , DR, QAZ , QZ , HEAT 

COMMON/ BLOKD/TGAS ( 50 ) , TO ( 15 , 10 ) , T ( 15 , 10 ) , AREA , GASMAS , 

, TGASO ( 50 ) , GASMIX , TGASIN , DP , VF , DPHI , CPGAS , 

, DENGAS , TP ( 15 , 10 ) ,TGASP(50) , TGAST ( 50 , 100 ) 

COMMON/ BLOKF/ I , J , L, NI , NJ , NL, ITIME , MODE 

COMMON/ BLOKG/KE (15,10) ,KW(15, 10) , KN , K3 , HWGAS ( 10 ) ,HWSOL 
COMMON/BLOKJ/A ( 15 ) , B ( 15 ) , C ( 15 ) , D ( 15 ) 

C 

INTEGER DTIME 
REAL KN , KS , KE , KW 
C 

CVPURE=RMELT ( I , J ) *CVLIQ+ ( 1 . 0-RMELT ( I , J) ) *CVSOL 

DENPURE=DENLIQ 

DR=RO/ (NI-1) 

RS= ( NI-I ) *DR 
C 

IF(I.EQ.l) THEN 
DR=DRM 

DRS= ( (RO/ (NI-1 ) )+DR)/2 
GO TO 100 

ELSEIF (I.EQ.2) THEN 
DRN= ( DR+DRM) /2 
DRS=DR 
GO TO 100 

ELSEIF (I.EQ.NI) THEN 
RS=0 . 0 
DRN=DR 
GO TO 100 
ELSE 

GO TO 90 
ENDIF 
C 


90 


DRN=DR 
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DRS=DR 

100 RN=RS+DR 

RE=RS+DR/2 

RW=RE 

DV= ( DPHI / 2 ) * ( RN-RS ) * ( RN+RS ) 

PF=1-VF 

AE=KE(I, J) *DR/ (RE*DPHI ) 

AW=KW(I, J) *DR/ (RW*DPHI ) 

C 

IF (I.EQ.NI) THEN 
C(I)=0.0 
GO TO 250 
ENDIF 
C 

AWALL= PI* RO / N J 

C(I) =- (KS*RS*DPHI/DRS ) 

C 

IF (I.EQ.l) GO TO 290 
C 

250 A ( I ) =- (KN*RN*DPHI/DRN) 

BC=H*SAPUV*DV*TGAS (L)+PF*DENPURE*CVPURE*DV*TO(I, J)/DTIME 
/+PF*DENPURE*HFL*DV* (ROMELT ( I , J) -RMELT(I, J) ) /DTIME 
B ( I ) =AE+AW-C ( I ) -A ( I ) +H*SAPUV*DV+PF*DENPURE*CVPURE*DV/ DTIME 
D ( I ) =AE*T ( I , J- 1 ) +AW*T ( I , J+ 1 ) +BC 
GO TO 300 
C 

290 C(I)=C(I) 

B(I) =AE+AW-C(I) +DENMET*CVMET*DV/DTIME+HWGAS ( J) *AWALL 
CALL TSOLAR 

BC= (DENMET*CVMET*DV*TO ( I , J) /DTIME) +QAZ+HWGAS ( J) *AWALL* 
*TGAS (L) 

D ( I ) =AE*T ( I , J-l ) +AW*T ( I , J+l) +BC 
A(I) =0. 0 
C 

300 RETURN 
END 




CIO 






c ********************************************************* 

Q Qle ic ie ic -k * ic *k * * * * * 1c * * * ic ie * * ic ic * * * * ic ic -k * it * * *k -k vi ******************** 

c c * 

C C SUBROUTINE CONDUCT * 

C C * 

C C This program calculates the interface conductivities to* 

C C the north, south, east, and west of the element (I,J). * 

C C * 

Q ( 2 ********************************************************* 

Q C********************************************************* 

C C 

SUBROUTINE TCON 
C 

COMMON/ BLOKF/ I , J , L, NI , NJ , NL, ITIME , MODE 

COMMON/ BLOKA/RMELT (15,10), KLIQ , KSOL, KGAS , KMET , VOID , EPS , 

, DRM , RO , DENLIQ , CVLIQ , CVSOL , DENMET , CVMET , HFL, 

, SAPUV , PI , PR , REP , H , CVPURE , DENPURE , DENSOL , 

, ROMELT (15,10) , RPMELT ( 15 , 10) 

COMMON/ BLOKD/TG AS (50) , TO ( 15 , 10 ) , T ( 15 , 10 ) , AREA, GASMAS , 

, TGASO ( 50 ) , GASMIX, TGASIN, DP, VF, DPHI , CPGAS , 

, DENGAS , TP (15,10) ,TGASP(50) , TGAST ( 50 , 100 ) 

COMMON/ BLOKG/KE (15,10) ,KW(15,10) , KN , KS , HWGAS ( 10 ) ,HWSOL 
COMMON/ BLOKO/KPSOL, KPURE , K ( 3 , 3 ) 

C 

REAL KPSOL , KPURE , DFT , DF , KSTAG , KCOND , KRAD , KRO , FS , KWKF 
REAL KFSTAG , KFF , KWALL , KWRAD , K , KLIQ , KSOL , KGAS , KMET , KWGRAD 
REAL KN , KS , KE , KW , NUM 
C 

DO 110 11=2,3 

IF ( (I.EQ.l) .AND. (II.NE.3) ) GO TO 110 
C 

IF ( (I.EQ.NI) .AND. (II.EQ.3) ) THEN 
KN=KS 
KS=0 . 0 
GO TO 155 
ENDIF 
C 

DO 120 JJ=2 , 3 
C 

IF (II .EQ. 2) THEN 

IF (JJ .EQ. 2) THEN 

VOID=(l. O-DENLIQ/DENSOL) * ( 1 . 0-RMELT ( I , J) ) 

ELSE 

VOID=(l. O-DENLIQ/DENSOL) * ( 1 . 0-RMELT ( I , J+l) ) 

ENDIF 

ELSE 

VOID= (1. O-DENLIQ/DENSOL) *(1. 0-RMELT (1+1, J) ) 

ENDIF 

C 

IJ=J 




Cll 


c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


IF (II.EQ.2) IJ=J-2+JJ 
IF ( (J.EQ.NJ) .AND. (JJ.EQ.3) ) THEN 
KW ( I , J) =0 . 0 
GO TO 105 
ENDIF 
C 

C********************************************************* 

C * 

C CALCULATE THE EFFECT OF VOID ON CONDUCTIVITY OF SOLID * 
C USE CORRELATION BY MEREDITH & TOBIAS TO CORRECT * 

C CONDUCTIVITY * 

C * 

C***************************■********'*•^t**'*■***************** 

c 

KPSOL=KSOL* (2 . 0* (1 . O-VOID) +0 . 6135*VOID**2 .333-1.6* 

*VOID**3 .333)/ 

/ (2 . O+VOID+O . 6135*VOID**2 . 333-0 . 6975*VOID**3 .333) 

C 

C********************************************************* 

C FIND CONDUCTIVITY OF PURE MATERIAL * 

c********************************************************* 
C 

JI=I 

IF (II .EQ. 3) JI=I+1 

KPURE=1 . 0/ ( (RMELT( JI , IJ) /KLIQ) + ( ( 1 . 0-RMELT ( JI , IJ) )/KPSOL) ) 

HRS=0. 227* (EPS/ (2 . 0-EPS) ) * (T ( I , J) **3) *10E-9 

HRV= . 2 2 7 *T ( I , J ) **3*10E-9/ (1+VF* (1. 0-EPS)/ (2* ( 1-VF) *EPS) ) 

C 

C* ******************************************************** 

c * 

C FIND THE EFFECTIVE CONDUCTIVITY OF THE BED USING METHOD* 

C * 

C OF SUMMED RESISTANCE OVER VARIOUS HEAT TRANSFER PATHS * 

C * 

C " HANDBOOK OF HEAT TRANSFER APPLICATIONS, * 

C 2nd EDITION, 1985, pp. 6-10,6-18. * 

Q-k**** ******************************************** ******** 

C * 

DFT- ( (KPURE-KGAS) /KPURE) **2 

DF=( (0.25* DFT ) / ( ALOG ( KPURE/ KGAS ) -DFT) ) -KGAS/ ( 3 * KPURE ) 

NUM— ( (KPURE-KGAS) /KPURE) **2 

Dl=ALOG( (KPURE/KGAS) -.9128* ( ( KPURE/KGAS ) -1.0)) 

D2= ( (KPURE-KGAS) /KPURE) * . 0872 

PHI=( (1.0/12) *NUM/ ( D1-D2 ) ) -(2.0/3) * (KGAS/KPURE) 

D= ( 1/PHI ) + ( DP*HRS/KGAS ) 

KSTAG=VF* ( 1+HRV*DP/ KGAS ) + (1-VF) / 

/( (1. 0/D) +(2.0/3) *KGAS/KPURE) 

KCOND= (KSTAG+0 . 11*PR*REP) *KGAS 
K(II , JJ) =KCOND 
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C 

C 

C 

C 

C 

C 

C 

C 

c 

c 


105 IF (II.EQ.3) GO TO 110 
120 CONTINUE 
110 CONTINUE 

C * 

C ***** * ********* * * * ************ * * * * * * * * ************* * ***** 

C FIND WALL-TO-BED CONDUCTIVITY USING THE METHOD OF * 

C * 

C "EFFECTIVE PATH LENGTH" * 

C * 

C "HANDBOOK OF HEAT TRANSFER APPLICATIONS 2nd EDITION, * 

C 1985, pp. 6-10,6-18. * 

q*** ****************************** *********** ************* 

c * 

IF (I.EQ.l) THEN 

190 KWKF=2 . 0*0 . 7+ ( 1 . 0-0 . 7) / (DF+ (KGAS/ ( 3 . 0*KPURE) ) ) 

KFSTAG=1 . 0/KWKF- ( 0 . 5*KGAS/KS TAG ) 

KFF=1 . O/KFSTAG+O . 054 *PR*REP 
KWALL=KFF*KGAS/DP 

DEN= ( VF* ( 1 . 0-EPS ) *10E+6) / (2 . 0* ( 1-VF) *EPS) 

KWGRAD= (0.227*T(I,J) **3/ (DEN+1 . 0) )/1000 
HWGAS (J) =KWALL 

HWSOL=( (0.227* (EPS/ (2. 0-EPS) ) *(T(I,J) **3*10E-6) )/1000)+ 
+ KWRAD 

KE ( I , J) =KMET 

IF (J .EQ. 1) KE ( I , J) =0 . 0 
KW(I, J)=KMET 

IF (J .EQ. NJ) KW ( I , J) =0 . 0 
FS=DRM/ (DRM+RO/ (NI-1 . 0) ) 

K(3 , 2) =KWKF*KGAS+0. ll*PR*REP*KGAS+fIWSOL* 

* (RO/ (NI-1 . 0) +DRM) /2 . 0 

KS=1 . 0/ ( FS/KMET+ ( 1 . 0-FS ) /K (3,2) ) 

KN=0 . 0 
GO TO 200 
ENDIF 


KN=KS 

KS=1 . 0/ ( 0 . 5/K ( 2 , 2 ) +0 . 5/K ( 3 , 2 ) ) 


155 IF (J.EQ.l) THEN 
KE (I , J) =0 . 0 
GO TO 195 
ELSE 

KE ( I , J) =KW ( I , J-l ) 
ENDIF 


IF ( (KW (I , J) . EQ. 0) .AND. (J .EQ. NJ) ) GO TO 200 
195 KW(I,J)=1.0/(0.5/K(2,2)+0.5/K(2,3) ) 

200 RETURN 
END 
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C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


0 ********************************************************** 
c * 

C SUBROUTINE GASPROP * 

C * 


C This program calculates the effective thermal properties* 
C and convective coefficients for the gas coolant. * 
C The gas coolant is assumed to be a mixture of He and Xe.* 
C Properties of the two constituents are included as curve* 
C fits with the corresponding coefficients included as * 
C data in the routine. * 


C 


* 






c 

C* ********************************************************* 

c * 

C REFERENCES * 

C M. R. Vanco, "Analytical Comparison of Relative Heat * 
C Transfer Coefficients and Pressure Drops of Inert * 

C Gasesand Their Binary Mixtures", Lewis Research Center* 


C * 
C A. J. Chapman "Fundamentals of Heat Transfer" * 
C Macmillan, 1978, p 329. * 
C * 
c********************************************************** 

c 

c********************************************************** 
c * 

C EQUATIONS * 
C K1 = Kll + K12* (TGAS - K13)**K14 * 
C K2 = K21 + K2 2 * ( TGAS - K23)**K24 * 
C * 
C VISC1 = VISC11 + VISC12* (TGAS - VISC13 ) **VISC14 * 
C VISC2 = VISC21 + VISC22* (TGAS - VISC23 ) **VISC24 * 
C * 
C **** ****************************************************** 


c 

SUBROUTINE TGAS2 (KGAS , H, PR, REP) 

REAL K, K1 , K2 , Kll , K12 , K13 , K14 , K21 , K22 , K23 , K24 , MW, MW1 , MW 2 , 
, KGAS 

COMMON/ BLOKD/ TGAS (50) , TO ( 15 , 10 ) , T ( 15 , 10 ) , AREA , GASMAS , 

, TGASO ( 50 ) , GASMIX, TGAS IN, DP, VF, DPHI,CPGAS, 

, DENGAS , TP (15,10) ,TGASP(50) , TGAST ( 50 , 100 ) 

COMMON/ BLOKF/ I, J, L,NI,NJ,NL, ITIME , MODE 

COMMON/ BLO KG/ KE ( 15 , 10) ,KW(15,10) , KN , KS , HWGAS ( 10 ) ,HWSOL 
DATA GAMMA, PI, RUNIV , MW1 , MW2/ 

/l. 66, 3. 14 159, 8. 3 14 34, 4. 0,13 1.0/ 

DATA Kll, K12,K13,K14/0. 1580E-3 , 0 . 9086E-6 , 333 .333, 0.794/ 


C14 


C 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 


DATA K2 1 , K22 , K2 3 , K24/0 . 64 0E-5 , 0 . 0454E-6 ,333.333,0.800/ 

DATA VISC11 , VISC12 , VISC13 , VISC14/ 

/ . 2158E-4 , .1022E-6, 333.333, .825/ 

DATA VISC21, VISC22 , VISC23 , VISC24/ 

/0.253E-4,0.2452E-6, 333. 333, 0.77/ 

FI (X, Y) = (1 . 0+2 .41* (X-Y) * (X-0 . 142*Y) / (X+Y) **2) 

F2 (X1,X2,X3,X4,X5,X6)= 

= ( (1.0+(X1/X2) **0.5* (X3/X4) **0.25) * * 2 ) / 

/(2. 0*2. 0**0. 5*(1. 0+X5/X6) **0.5) 

C 

0 ********************************************************** 

c * 

C FIND THE THERMAL CONDUCTIVITY & VISCOSITY OF PURE GASES * 
C * 

0 ********************************************************** 

c 

K1=K11+K12 * (TGAS ( L) -K13 ) **K14 
K2=K21+K22* (TGAS (L) -K23) **K24 
VISC1=VISC11+VISC12 * (TGAS (L) -VISC13) **VISC14 
VISC2=VISC21+VISC22* (TGAS (L) -VISC23) **VISC24 
C 

0 ********************************************************** 
C * 

C FIND THE WEIGHTING PARAMETERS FROM REFERENCE * 

C * 

c 

PHI12=F2 ( VISC1 , VISC2 , MW 2 , MW1 , MW1 , MW2 ) 

PHI 2 1=PHI 12* (VISC2*MW1) / (VISC1*MW2) 

PSI12=F2 ( K1 , K2 , MW1 , MW2 , MW1 , MW2 ) *F1 (MW1 , MW2 ) 

PSI2 1=F2 (K2 , K1 , MW2 , MW1 , MW2 , MW1 ) *F1 (MW2 , MW1 ) 

C 

q* ************************************ ******************** * 
C * 

C SOLVE FOR THE THERMAL CONDUCTIVITY & VISCOSITY OF MIX * 
C * 

0 ********************************************************** 

c 

VISC=VISC1/ ( 1 . 0+PHI12* ( 1 . 0-GASMIX) /GASMIX) + 

/ VISC2/(1.0+PHI21* (GASMIX/ (1. 0-GASMIX) ) ) 

KGAS=K1/ ( 1 . 0+PSI12* ( 1 . 0-GASMIX) /GASMIX) + 

/ K2/ (1.0+PSI21* (GASMIX/ (1. 0-GASMIX) ) ) 

0 ************************************* ******************** 

c * 

C SOLVE FOR SPECIFIC HEAT OF GASMIX * 

C * 

c 


MW=GASMIX*MW1+ ( 1 . 0-GASMIX) *MW2 
CPGAS=GAMMA*RUNIV/ ( (GAMMA-1.0) *MW) 
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C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 


c 

c** ************************************************* ******* 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


* 

CORRELATION RECOMMENDED FOR GAS FLOW IN A BED OF SPHERES* 


_0 . 575 * 

(ej) = (ej) = 2.06 (Re) * 

h m * 

* 

j are the Colburn j factors, and e is the void fraction * 

* 

REFERENCE * 

Incropera, Frank P. , "Fundamentals of Heat and Mass * 

Transfer", 2nd Edition, John Wiley & Sons Inc., 1981, * 

p 292,352. * 

* 


C************************** ************************** ****** 


C 


PR=CPGAS*VISC/ (KGAS) 

REP=GASMAS*DP/ (VISC*AREA) 

H= ( 2 . 06/VF) * (KGAS/DP) *PR** ( 1 . 0/3 . 0) *REP**0 . 425 
DENGAS=11. 22 
C 

RETURN 

END 
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C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


c******************************* *************************** 
c********************************************************** 

c * 

C SUBROUTINE GASTEMP * 

C * 

C This program solves for the temperature of the working * 
C fluid. The energy equation has been discretized using * 

C the upwind scheme to calculate the convection term. * 

C Ref: " Numerical Heat Transfer And Fluid Flow ", * 

C Suhas V. Patankar. * 

C * 

Q* ******* ******* ******* ******* ***************************** 

q********************************************************** 

c 

SUBROUTINE TGASTEMP 

COMMON/ BLOKA/RMELT (15,10), KLIQ , KSOL, KGAS , KMET , VOID, EPS , 

, DRM , RO , DENLIQ , CVLIQ , CVSOL , DENMET , CVMET , HFL, 

, SAPUV , PI , PR , REP , H , CVPURE , DENPURE , DENSOL , 

, ROMELT (15,10) , RPMELT (15,10) 

COMMON/ BLOKD/TGAS (50) , TO (15, 10) ,T(15, 10) , AREA , GASMAS , 

, TGASO ( 50 ) , GASMIX, TGASIN , DP, VF , DPHI , CPGAS , 

, DENGAS , TP ( 15 , 10 ) ,TGASP(50) , TGAST ( 50 , 100 ) 

COMMON/ BLOKB/DZ , DTIME , DR, QAZ , QZ , HEAT 
COMMON/ BLOKF/ I , J , L, NI , NJ , NL, ITIME , MODE 

COMMON/ BLOKG/KE ( 15 , 10 ) ,KW(15,10) , KN , KS , HWGAS ( 10 ) ,HWSOL 

REAL KGAS 
INTEGER DTIME 


DR=RO/ (NI-1) 

QWALL1=0 . 0 
QWALL2=0 . 0 
QBED1=0. 0 
QBED2=0 . 0 

AWALL=(PI*RO/NJ) *DZ 
DO 11 JQ=1 , NJ 
DO 11 IQ=2 , NI 

RS= (NI-IQ) *DR 
RN=RS+DR 

DV= (DPHI/2 ) * (RN-RS) * (RN+RS) *DZ 
IF (IQ .EQ. 2) THEN 

QWALL1=QWALL1+2*HWGAS (JQ) * AWALL*T ( 1 , JQ ) 
QWALL2=QWALL2+2*HWGAS (JQ) *AWALL 
ENDIF 


QBED1=QBED1+2*H*SAPUV*DV*T(IQ, JQ) 
QBED2=QBED2+2*H*SAPUV*DV 
11 CONTINUE 


VOLUME= ( PI*R0**2 ) *DZ 
AB=GASMAS*CPGAS 

B= (QWALL1+QBED1) +DENGAS*VF*CPGAS*TGASO (L) *VOLUME/DTIME 
A=AB + (DENGAS*VF*CPGAS*VOLUME/DTIME) +QWALL2+QBED2 

IF (L.EQ.l) THEN 
TB=TGASIN 

ELSE 

TB=TGAST ( L-l , ITIME ) 

ENDIF 

TGAS (L)=(AB*TB+B)/A 

RETURN 

END 
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C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


c 


c 


c 


c 


c 


Q********* ***************************** ************** ****** 

c * 

C SUBROUTINE MELT * 

C * 

C * 

C * 

c * 

c * 

c 


SUBROUTINE TMELT ( THMELT ) 

COMMON/ BLOKF/ I , J , L , NI , NJ , NL , ITIME , MODE 

COMMON/ BLOKA/RMELT (15,10), KLIQ, KSOL, KGAS , KMET , VOID, EPS , 

, DRM , RO , DENLIQ , CVLIQ , CVSOL, DENMET , CVMET , HFL , 

, SAPUV , PI , PR , REP , H , CVPURE , DENPURE , DENSOL , 

, ROMELT ( 15 , 10) , RPMELT ( 15 , 10) 

COMMON/ BLOKD/TGAS (50) ,TO(15, 10) ,T(15, 10) , AREA , GASMAS , 

, TGASO ( 50 ) , GASMIX , TGASIN , DP , VF , DPHI , CPGAS , 

, DENGAS , TP (15, 10) ,TGASP(50) , TGAST ( 50 , 100 ) 

COMMON/ BLOKB/DZ , DTIME , DR, QAZ , QZ , HEAT 
COMMON/ BLOKX/FCTR 

DO 10 11=1, NI 

T(II, J)=T( II, J) -THMELT 
10 CONTINUE 

DO 11 1=2, NI 
TN=T(I-1, J) 

IF (I .EQ. NI) THEN 
TS= (T ( I , J ) ) 

ELSE 

TS=T ( 1+1 , J) 

ENDIF 

RMELT ( 1 , J) =1 . 0 

DMELT=FCTR*T ( I , J ) *CVPURE/HFL 

IF ( (DMELT . LT. 0.0) .AND. (QAZ .LT. 0.0)) THEN 

IF (T ( I , J) .LE. 0.0) RMELT (I, J) =RMELT ( I , J) +DMELT 
ELSEIF (I .EQ. 2) THEN 

IF ( (T ( 1 , J) .GE. 0.0) .AND. (QAZ .GE. 0.0)) THEN 
RMELT ( I , J) = RMELT ( I , J) + DMELT 
ENDIF 

ELSEIF (T ( I , J) .GE. 0.0) THEN 





RMELT ( I , J ) =RMELT ( I , J ) +DMELT 
ELSE 

RMELT ( I, J) =0.0 
ENDIF 

RMELT ( I , J) =AMAX1 (0.0, RMELT (I , J) ) 
RMELT ( I , J ) =AMIN1 (1.0, RMELT ( I , J ) ) 
11 CONTINUE 

DO 13 1=1, NI 

T ( I , J) =T (I , J) +THMELT 
13 CONTINUE 

RETURN 

END 
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I 

f 

t 


i 


i 


» 


I 


I 


I 


I 


c************************************************************** 


c * 

C SUBROUTINE SOLAR * 
C * 
C This program calculates the solar incidence on * 
C the outside perimeter of the containment ring * 
C * 
C q" = q" (direct) *cos (angle) + q" (indirect) * 
C * 


c*************************************************************** 

c*************************************************************** 

c 

SUBROUTINE TSOLAR 
C 

DIMENSION CS (21) 

REAL LENGTH 

COMMON/ BLOKF/ I , J , L, NI , NJ , NL, ITIME , MODt 

COMMON/ BLOKB/DZ , DTIME , DR, QAZ , QZ , HEAT 

COMMON/ BLOKS/Q1 DIR, Q1IN, Q2DIR, Q2IN 

DATA CS(1) ,CS(3> ,CS(5) ,CS(7) ,CS(9) ,CS(11) ,CS(13)/ 

/l. 0,6*2. 0/ 

DATA CS (15) , CS ( 17 ) ,CS(19) ,CS(2) ,CS(4) ,CS(6) ,CS(8)/ 

/3 *2 .0,4*4 . 0/ 

DATA CS(10) , CS ( 12 ) ,CS(14) ,CS(16) ,CS(18) ,CS(20) ,CS(21)/ 
/6*4. ,1./ 

DATA PI/3.14159/ 

DATA QZ1,QZ2,QZ3, QZ4/24 .28,0.08,0.75,0.26244/ 

C 

Q****************************************** ********************* 


c * 

C QZ USES SIMPSON ' S 1/3 INTEGRATION TO FIND THE * 
C INTETRATED AVERAGE VALUE OF THE AXIAL POWER CURVE * 
C OVER THE INTERVAL DZ * 
C * 


0 *************************************************************** 

c 

LENGTH=0. 91-4 4 
POS=DZ*L 

XL= (POS-DZ ) /LENGTH 
QZ=0 . 0 
c 

DO 15 11=1,21 

IF (II .GT. 1) XL=XL+0 . 05*DZ/LENGTH 
QZZ=QZ1*XL*EXP( 1 . 0- (XL/QZ2) **QZ3) +QZ4* (XL) **0.5 
QZ=QZ+CS (II) * (QZZ)/60. 0 
15 CONTINUE 
C 

q * ********* ****************’****** ****************************** * 

c * 


I 


o o o o 


c 21 


FIND THE PROJECTION OF THE DIRECT SOLAR INCIDENCE * 

* 

*************************************************************** 


ANGLE 1=PI* FLO AT (J) /FLOAT (NJ) 

ANGLE2=PI*FLOAT (J-l) /FLOAT (NJ) 

C 

IF (J+J-2 .GT. NJ) THEN 
SINANGLE=0 . 0 

ELSEIF (J+J .GT. NJ) THEN 
SINANGLE=1 . O-SIN (ANGLE2 ) 

ELSE 

SINANGLE=SIN ( ANGLE 1) -SIN ( ANGLE2 ) 

ENDIF 

C 

c** **************************************************** ********* 
C * 

C ADD DIRECT AND INDIRECT INCIDENT SOLAR RADIATION * 

C * 

C***********************************************************'**** 

c 

IF (MODE .EQ. 1) THEN 
QDIR=Q1DIR 
QINDIR=Q1IN 
ELSE 

QDIR=Q2DIR 

QINDIR=Q2IN 

ENDIF 

DPHI=PI/NJ 

ROCAN=0. 04191+0. 00127 

QAZ= (QDIR*SINANGLE + QINDIR*DPHI) *ROCAM*QZ 
IF ( J .EQ. 1 ) HEAT=0 . 0 
HE AT=HE AT+ 2 . 0*QAZ 
C 

RETURN 

END 
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C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 


c * ********************************************************* 
c * ********************************************************* 
C SUBROUTINE SOLVER * 

C This program solves a tridiagonal matrix using an * 

C elimination method * 

C * 

C [ A] T ( 1-1 ) + [B]T(I) + [C]T(I+1) = [D] * 

C * 

£*★******★★*****★★★■*****★*******★***★* .'. .h-k'k'k'k'kick'k'k'k'kie-kicrkic-k'k-k 

c********************************************************** 

c 

SUBROUTINE TSOLVER 


COMMON/ BLOKD/TGAS (50) ,TO(15, 10) ,T(15, 10) , AREA , GASMAS 
, TGASO (50) , GASMIX, TGASIN, DP, VF , DPHI , CPGAS , 

, DENGAS , TP (15,10) ,TGASP(50) ,TGAST(50, 100) 

COMMON/ BLOKJ/ A ( 1 5 ) ,B(15) ,C(15) ,D(15) 

COMMON/ BLOKF/ I , J , L, NI , NJ , NL, ITIME , MODE 
C 

Q********************************************************** 


c * 

C PROCEED WITH GAUSSIAN ELIMINATION TO DEVELOP * 
C AN UPPER DIAGONAL MATRIX * 
C * 


C 

DO 1 1=1 , NI-1 

R=A ( 1+1 ) /B ( I ) 

B ( 1+1) =B ( 1+1) -R*C(I) 

D ( 1+1 ) =D ( 1+1 ) -R*D ( I ) 

1 CONTINUE 
C 

C ************************************** ******************** 

c * 

C BEGIN BACK SUBSTITUTION * 

C * 

C 

T(NI,J)=D(NI)/B(NI) 

DO 2 I=NI-1 , 1,-1 

T(I,J)=(D(I)-C(I)*T(I+1,J) )/B(I) 

2 CONTINUE 
RETURN 
END 
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program main 

c 

common/bloka/a (15) , b ( 15) , c ( 15) ', d ( 15) 
common/blokb/ni , ne ( 10) , j , t ( 15 , 10) 

common/blokc/thmelt , tgasin , gasmas , cpgas , ntu , deltp , 

+ deltdz 

common/blokd/tgasp , gasmix , ro ( 15 ) , ht , w , dxx , p 
common/bloke/ rmelt (15,10) , deni iq, densol , walls , cvmet , 

+ denmet , kmet , kl iq , ksol , cvl iq , cvsol 

common/blokf/dz,dtime,rpo(10) , qaz , hfl , romelt ( 15', 10) ,cl,wf , 
+ to (15, 10) , rpmelt (15, 10) , dx, dy , th, qt , le (15) , 

+ apoo ( 15 , 10) , apon ( 15 , 10) 

common/blokg/k, kn, ks, ke, kw,den, cv 
dimension cs(21) ,q(20) ,tp(15,10) , thgas (20,20) , 

+ theta (20, 10) ,thn(15, 10) , rmel (15, 10) , nn(10) 

real ntu, length, kliq, ksol , kmet , k, kn , ks , ke , kw, den , cv, le 
c 

data cs(l) ,cs(3) ,cs(5) ,cs(7) ,cs(9) ,cs(ll) ,cs(13)/ 

/l. 0, 6*2 . 0/ 

data cs ( 15) ,cs(17) ,cs(19) ,cs(2) ,cs(4) ,cs(6) ,cs(8)/ 

/3*2.0, 4*4.0/ 

data cs(10) ,cs(12) ,cs(14) ,cs(16) ,cs(18) ,cs(20) ,cs(21)/ 

/ 6*4 . , 1 . / 

open ( 8 , f ile= 1 g input . txt ’ ) 
open (7, file=' goutput.txt ' ) 
c 

read (8, *) imax, norbit , nj ,nl,ntime 

read (8, *) length, ht,w,th, si 

read (8 , *) gasmas , borl , cl , wf , tr 

read (8 , *) cvl iq, deni iq, kliq, walls, thme It 

read ( 8 , * ) cvsol , densol , ksol , errmax , hfl 

read ( 8 , * ) cvmet , denmet , kmet , gasmix , timer 

read (8 , *) qzl , qz2 , qz3 , qz4 , qz5 , qz6 

read ( 8 , * ) qd , qn , tgas , cor , ft 

p = 3.14 
ni = nj 
tgasp = tgas 
dy = ht/nj 

dz = length/float (nl) 
dtime = 60 . 0*timer/f loat (ntime) 
qtime = 54 . 0*60 . 0/dtime 
c 

do 2 j = 1 , nj 
ro(j) = ni*dy 
ne(j) = ni-j+2 
2 continue 
c 

do 100 1 = 1 , nl 

w = p* ( (nl-l+. 5) *dz/sl+tr)/cor 
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dx = w/nj 

dxx = (dx*dx+dy*dy) **0 . 5 
c 

do 10 jj = l,nj 
do 10 ii =1,15 
c 

if (ii .gt. ne(jj)) then 
t(ii,jj) = 0.0 
else 

t(ii,jj) = -0.00001 
endif 
c 

tp(ii,jj) = t(ii, j j ) 
rmelt (ii, j j ) = 0.0 
apon(ii,jj) = 0.0 
10 continue 
pos = l*dz 
xl = pos/ length 

qzz = qzl+qz2*xl+qz3*xl*xl+qz4* (xl) **3 . 0+qz5* (xl) **4 . 0 
+ +qz6* (xl) **5. 0 

qz = qzz-qzzp 
qzzp = qzz 
c 

do 100 no = l,norbit 
do 100 n = l,ntime 
c 

if (1 .eq. l) thgas(n,no) = tgas 
c 

do 30 jj = l,nj 
do 30 ii = 1,15 
c 

if (ii.gt.ne( j j ) ) then 
t ( ii , j j ) =o.o 
apon(ii,jj) = 0.0 
endif 
c 

to(ii, j j) = t ( ii , j j ) 
romelt(ii, j j) = rmelt ( ii , j j ) 
apoo(ii, j j) = apon ( ii , j j ) 

30 continue 
c 

if (n.gt . qtime) then 
q(n) = qn 
else 

q(n) = qd 
endif 
c 

qaz = q(n)*qz/nj 
icount = 0 
35 iflag = 0 



icount = icount+1 
tgasin = thgas(n,no) 
ta=0 . 0 
tt=0 . 0 

do 24 j = l,nj 
tt = tt + ne( j) -1 
le ( j ) = ro(j)-(j-l) *dy 
f = ( 1-exp ( -cl* le(j ) ) ) 
ta = ta+f 
24 continue 

qt = qaz* (nj -ta) * ( 1-wf ) 

qt = qt/tt 

qlat = 0.0 

qsens = 0.0 

qrad = 0.0 

tac = 0 

heat = 0.0 

hd = 0.0 

do 80 j = l,nj 
rpo(j) = ro ( j ) 

80 continue 

do 40 j = l,nj 
call gasprop (dz , nj ,h, 1) 
heat = heat + qaz*dz 
hd = hd+h 

call conduct (borl , i , n j , icount , ft , qlat , qsens , qrad , tac) 
call heatgas ( n j ) 
call solver(i) 

40 continue 

hd = hd/nj 

do 50 jj = l,nj 
do 50 ii = l,ne(jj) 
error = abs (t ( ii, j j ) -tp ( ii, j j ) ) 
if (error ,gt. errmax) if lag = 2 
tp(ii,jj) = t ( ii , j j ) 
rpraelt(ii, j j ) = rmelt ( ii , j j ) 

50 continue 

qgas = gasmas*cpgas*deltdz 
qtotal = qlat+qsens+qgas 
qerror = abs ( 1 . 0-qtotal/heat) 

if (n .gt. qtime) then 

if (qerror .gt. 0.01) if lag = 


2 


D4 
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else 

if (qerror .gt. 0.001) iflag = 2 
endif 

tgasp = thgas (n, no) +deltp 

if (iflag .gt. 1) then 

if (icount ,ge. imax) then 
write (7,1000) icount 
else 

goto 35 
endif 
endif 


thgas (n, no) = thgas (n, no) +deltdz 

do 47 j =l,nj 
47 continue 

if (no .It. norbit) goto 100 
write (7 , 1050) 
write (7 , 1250) 
write (7 , 1100) 1 , no, n 
write (7 , 1150) heat , qz 
write (7 , 1151) qtotal 
1151 format (lx, 'qtotal = ' ,flo.8) 
nr = 0.0 

do 83 j = l,nj 
nn( j )=ne(j )+j-l 
83 continue 

do 34 j = i,nj 
nr = amaxl(nr,nn( j) ) 

34 continue 

do 36 j = l,nj 
do 37 i = l,nr 

if (i.gt.ne(j)) then 
rmelt(i,j) = 0.0 
theta (i,j) = 0.0 
else 

theta (i,j) = t (i, j ) +thmelt 
endif 

37 continue 
36 continue 


do 85 j = l,nj 
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do 85 i = l,nr 
c 

if (j.ne.l) then 

thn( (i+j-l) , j ) = theta (i,j) 
rmel ( (itj -1) , j ) = rmelt(i,j) 
else 

thn(i,j) = theta (i,j) 
rmel(i,j) = rmelt(i,j) 
endif 

c 

85 continue 
c 

do 70 i = nr, 1,-1 

write (7 , 1200) i , (thn ( i, j ) , j=l, nj ) 

70 continue 
c 

write (7,1300) tgasp, thgas (n, no) , hd 
write (7, 1050) 
write (7, 1250) 
c 

do 82 i = nr, 1,-1 

write (7 , 1350) i , ( rmel (i,j) ,j=l,nj) 

82 continue 
c 

100 continue 

c************* ******************************************** ****** 

c * 

c input and output format statements * 

c * 

1000 format ( lx, ' failed to converge after ', i6 ,' iterations ' ) 

1050 format (lx,' ') 

1250 format ( lx , ' =============================================== 

===',' ==================== • ) 

1100 format ( lx, ' axial position = ',il,' orbit = ',il,' orbit 
& position = ' , i2 ) 

1150 format (lx, 'net solar heat input = ' ,lfl0.8,'qz = ' ,lf9.4) 
1200 format (lx, ' i = ',i2,20f7.1) 

1300 format ( lx, ' avg gas temp = ' ,f6.1,'k out gas temp ', 
,f6.1,'kh = 1 , f7 . 4 , 'kw/m2-k' ) 

1350 format (lx, ' i = ' ,i2,20f7.4) 
close(8) 
close(7) 
stop 
end 


subroutine gasprop (dz , nj , h, 1) 

real k, kl , k2 , kll , kl2 , kl3 , kl4 , k21 , k22 , k23 , k24 , kgas , 

+ mw,mwl,mw2 , ntu,nu 
common/blokb/ni , ne (10) , j , t (15 , 10) 

common/blokc/thmelt , tgasin , gasmas , cpgas , ntu , deltp , 

+ deltdz 

common/blokd/ tgasp, gasmix, ro( 15) , ht, w,dxx,p 

data gamma, pi , runiv,mwl , mw2/l . 66,3.14159,8.31434,4,131/ 
data kll , kl2 , kl3 , kl4/0 . 1580e-3 , 0 . 9086e-6 ,333.333,0.794/ 
data k21,k22,k23,k24/0.6404e-5,0.0454e-6, 333 .333, 0.800/ 
data viscll, viscl2 , viscl3 , viscl4/ 

/ . 2 158e-4 , . 1022e-6 ,333.3, .825/ 
data visc21 , visc22 , visc23 , visc24/ 

/ . 253e-4 , . 2 4 52 e- 6 ,333.333, .77/ 
f 1 (x, y) = (1+2.41* (x-y) * (x-0 . 142*y) / (x+y) **2.0) 
f 2 (xl , x2 , x3 , x4 , x5 , x6 ) = 

= ( (1+(X1/X2) **.5*(X3/X4) **.25) **2.0)/ 

+ (2.0*2.0**.5*(l+x5/x6) **0.5) 

kl = kll+kl2*(tgasp-kl3) **kl4 
k2 = k21+k22* (tgasp-k23) **k24 

viscl = viscll+viscl2* (tgasp-viscl3) **viscl4 
visc2 = visc21+visc22* (tgasp-visc23) **visc24 
viscwl = viscll+viscl2* (t ( 1 , j ) +thmelt-viscl3) **viscl4 
viscw2 = visc21+visc22* (t ( 1 , j ) +thmelt-visc23) **visc24 

phil2 = f 2 (viscl , visc2 , mw2 , mwl , mwl , mw2 ) 
phi21 = phil2* (visc2*mwl)/ (viscl*mw2) 
psil2 = f 2 (kl , k2 , mwl , mw2 , mwl , mw2 ) *f 1 (mwl , mw2 ) 
psi21 = f 2 (k2 , kl , mw2 , mwl , mw2 , mwl) * f 1 (mw2 , mwl) 

vise = viscl/ (l+phil2* ( 1 - gasmix) /gasmix) + 

+ visc2/ ( l+phi21* (gasmix/ (1 - gasmix))) 

visew = viscwl/ (l+phil2* (1 - gasmix) /gasmix) + 

+ viscw2/ ( l+phi2 1* (gasmix/ ( 1 - gasmix))) 

kgas = kl/ ( l+psil2* ( 1 - gasmix) /gasmix) + 

+ k2/(l +psi21* (gasmix/ ( 1 - gasmix))) 

mw = gasmix*mwl+ ( 1-gasmix) *mw2 
cpgas = gamma*runiv/ ( (gamma-1) *mw) 

pr = cpgas*visc/ (kgas) 
dt = 2 . 0*ht*w/ (dxx*nj ) 
re = 2 . 0*gasmas/ ( (p* . 5*dt) *visc) 

nut = . 027*re** . 8*pr** . 333333* (visc/viscw) **. 14 

if ( (1 . eq. 1) . and. (re. It. 2100.0)) then 
nu = . 44* (dz/ (re*pr*dt) )**(-. 66) 


elseif (re .It. 2100.0) then 
nu = 2.35 

elseif ( (re.gt. 2100 . 0) . and. (re.lt. 10000 . 0 
nu = amaxl(2. 35, nut) 
else 

nu = nut 
end if 


h = kgas*nu/dt 


ntu = h*dxx*dz*nj/ 
if (l.eg.l) ntu = 
if (l.eq.2) ntu = 
if (l.eq.3) ntu = 
if (l.eq.4) ntu = 
if (l.eq.5) ntu = 


(gasmas*cpgas) 

.2337841 

.5330746 

.4739715 

.4153427 

.41995574 


) then 


return 

end 
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subroutine conduct 

+ (borl , i , nj , icount , ft , qlat , qsens , qrad , tac) 
c 

common/bloka/a ( 15) ,b(15) ,c(15) ,d(15) 
common/blokb/ni , ne (10) , j , t ( 15 , 10) 

common/blokc/thmelt , tgasin , gasmas , cpgas , ntu , deltp , 

+ deltdz 

common/blokd/ tgasp , gasmix , ro ( 15 ) , ht , w , dxx , p 
common/bloke/rmelt (15,10) , denliq, densol , walls , cvxnet , 

+ denmet , kmet, kliq, ksol , cvliq, cvsol 

common/blokf/dz , dtirae, rpo ( 10) , qaz,hfl,romelt(15, 10) ,cl,wf, 
+ to (15, 10) , rpmelt (15,10) , dx,dy , th, qt , le ( 15) , 

+ apoo ( 15 , 10) , apon ( 15 , 10) 

common/blokg/k, kn , ks , ke, kw , den, cv 
real k , kn , ks , ke , kw , kmet , ntu , ksol , kl iq , le 
c 

ne(j) = 200 
i = 1 
as = 0 
c 

1 if (j .eq. 1) then 
aw = 0 
tw = 0 

elseif (i.eq.l) then 
tw = t (i, j-1) 

else 

tw = t(i+l,j-l) 

end if 

c 

if (j .eq. nj) then 
ae = 0.0 
te = 0.0 

elseif (i.eq.l) then 
te = t (i , j+1) 
else 

te = t(i-l, j+1) 
endif 
c 

if (i .ge. 2) then 

r = ( float ( i) +j -2 . 5) *dy 

rn = r + .5*dy 
rs = r - .5*dy 

endif 

c 

c * ************************************************************** 

c set nodal constants for the wall * 

c * ************************************************************** 

c 

if (i .eq. 1) then 
an = kmet*dxx/th 
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if (j .ne. 1) aw = kmet*th/dxx 

if (j .ne. n j ) ae = kmet*th/dxx 

apon(i,j) = denmet*cvmet*th*dxx/dtime 

apo = apon(i,j) 

dapo =0.0 

dtapo = 0.0 

tcij = tgasin-thmelt 

hij = gasmas*cpgas* (l-exp(-ntu) )/ (dz*nj ) 
abw = wf* (exp (-cl*le ( j ) ) ) 
qr = qaz*abw 
qij = qr 
goto 40 
endif 
c 

Q* ******* ******* ******* ******* ********************************** 

c set nodal constants for outermost elements * 

q* ************************************************************* * 

c 

if (rpo(j) .gt. rn) then 

if (rpo(j) .le. (rn+dy) ) ne(j) = i+1 
endif 
c 

if ( ( j .eq. nj).and.(i . eq. 2)) then 
if (rpo ( j ) . le. rn) ne(j) = i 
endif 
c 

if (i .gt. (ni-j+2)) then 
rmelt(i,j) = 1.0 
c 

if (i .eq. ne(j)) then 
an = 0.0 

elseif (i .eq. (ne(j)-l)) then 

an = kliq*dx/ ( . 5* (rpo ( j ) -rn) + . 5*dy ) 
else 

an = kliq*dx/dy 
endif 
c 

if (i .eq. (ni-j+3)) then 

call property (cv, den, i-1 , j , ks) 
if (i.eq.ne(j)) then 

as = kliq*ks*dx/ (kliq* . 5*dy+ks* (rpo ( j ) -rs) * . 5) 
else 

as = kliq*ks*dx/ ( (kliq+ks) * . 5*dy) 
endif 

elseif (i .eq. ne(j)) then 

as = kliq*dx/ ( . 5* (rpo( j ) -rs+dy) ) 
else 

as = kliq*dx/dy 
endif 


c 



if (j .ne. nj) then 

if (rpo(j+l) .ge.rn) then 
if (i.eq. ne(j)) then 

ae = kliq*( (rpo( j)+rn) *.5-rs)/dx 
else 

ae = (kliq*dy)/dx 
endif 

elseif (rpo ( j+1) . le. rs) then 
ae = 0.0 

elseif (i.eq.ne(j)) then 

ae = (kliq* ( (rpo(j ) +rpo (j+1) )*. 5-rs) )/dx 
else 

ae = kliq* ( (rn+rpo( j+1) )*. 5-rs)/dx 
endif 
endif 

if (j .ne.l) then 

if (rpo ( j-1) .ge. rn) then 
i ; f (i.eq. ne(j)) then 

aw = kliq* ( (rpo (j ) +rn) *. 5-rs) /dx 
else 

aw = (kliq*dy)/dx 
endif 

elseif (rpo (j-1) . le. rs) then 
aw = 0.0 

elseif (i.eq.ne(j)) then 

aw = (kliq* ( (rpo (j ) +rpo (j-1) )*. 5-rs) ) /dx 
else 

aw = kliq*( (rn+rpo(j-l) ) *.5-rs)/dx 
endif 
endif 

if (i .eq. ne(j)) then 

apon ( i , j ) = denliq*cvliq* (1.0 - walls) * (rpo (j ) -rs) 
dx/dtime 

else 

apon(i,j) = denliq*cvl iq* ( 1 . 0 - walls) *dy*dx/dtime 
endif 

if (apoo ( i , j ) . eq. 0 . 0) apoo(i,j) =apon(i,j) 
if (qaz .gt. 0.0) then 
apo = apoo(i,j) 
dapo = apon(i,j) - apoo(i,j) 
dtapo = 0.0 
else 

apo = apon(i,j) 
dapo = 0.0 
dtapo = apoo(i,j) 
endif 

tcij = 0.0 


- apon ( i , j ) 
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hij = 0.0 

if (i .eq. ne(j)) then 

abm = 1 - exp (-cl* (rpo ( j ) -rs) ) 
else 

abm = (-exp (-cl* (rpo( j ) -rs) ) +exp(-cl* (rpo( j ) -rn) ) ) 
endif 

qr = qaz*abm + qt 
qij = qr 
goto 40 
endif 
c 

c** **************** ************************************* ******** 
c set nodal constants for the inner pcm * 

q* -k k * -k -k -k * * it * ic * * -k * * ic -k -kit 1c * * ic * *k -k -k 1e * -k -k * ic ic -k ic -k -k ic * * * * * * -k * -k -k * -k * * * -k * * -k k -k k 

23 call property (cv,den, i, j ,k) 

c 

if (i .eq. 2) then 

apon(i,j) = den*cv* . 5*dx*dy/dtime 
pcmmass = . 5*dx*dy*densol* ( 1-walls^ 
else 

apon(i,j) = den*cv*dx*dy/dtime 
pcmmass = dx*dy*densol* ( 1-walls) 
endif 
c 

if (apoo(i,j) .eq. 0.0) apoo(i,j) = apon(i,j) 
if (qaz.gt.0.0) then 
apo = apon(i,j) 
dapo = 0.0 

dtapo = apoo(i,j) - apon(i,j) 
else 

apo = apoo(i,j) 
dapo = apon(i,j) - apoo(i,j) 
dtapo = 0.0 
endif 
c 

tcij = 0.0 
hij = 0.0 

dmelt = t ( i , j ) *cv*borl/hf 1 
c 

if (i.eq. ne(j)) then 
an = 0.0 

elseif ( i . eq. (ni-j+2 ) ) then 
if ( ( i+1) . eq. ne ( j ) ) then 

an = k*kliq*dx/ (k* (rpo ( j ) -rn) * . 5+kliq* . 5*dy) 
else 

an = k*kliq*dx/ ( (k+kliq) * . 5*dy) 
endif 

else 

call property (cv, den, i+1 , j , kn) 
an = k*kn*dx/ ( (k+kn) * . 5*dy) 



D12 



endif 

if (i .eq. 2) then 

as = kmet*dxx/th 

else 

call property (cv, den, i-l,j , ks) 
as = k*ks*dx/ ( (k+ks) * . 5*dy) 

endif 

if ( j . ne . nj ) then 
if (i.eq.2) then 
ae = 0.0 
else 

call property (cv, den, i-1 , j+i , ke) 
ae = k*ke*dy/ ( (k+ke) * . 5*dx) 
endif 
endif 

•if (j.ne.l) then 

call property (cv, den, i+1, j-l,kw) 
aw = k*kw*dy/ ( (k+kw) * . 5*dx) 

endif 


c 

if (i.eq.2) then 

if (i.eq. ne(j)) then 

if (dmelt . gt . 0 . 0) then 

if ( (t(i-l, j) .ge.0.0) .or. (rmelt ( i+1 , j ) .ge. 1.0) 
& .or. (qaz.gt.0.0) ) 

& rmelt(i,j) = rmelt ( i , j ) tdmelt 

else 

if ( (t ( i-1, j ) .le.0.0) .or. (qaz.lt. 0.0)) 

& rmelt (i,j) = rmelt ( i , j ) +dmelt 

endif 


else 

if (dmelt .gt. 0 . 0) then 

if (rmelt(i+l, j) .ge. 1.0) 

& rmelt(i,j) = rmelt ( i , j ) +dmelt 

else 

if ( (t(i-l, j) .le.0.0) .or. (qaz.lt. 0.0)) 

& rmelt (i,j) = rmelt ( i , j ) +dmelt 

endif 
endif 

elseif (i.eq.ne(j)) then 
if (dmelt .gt. 0.0) then 

if (qaz .gt. 0.0) rmelt(i,j) = rmelt ( i , j ) +dmelt 
else 



if ( (rmelt(i-l, j ) .le. 0.0). or. (qaz .It. 0.0)) 
rmelt(i,j) = rmelt ( i ,j ) +dmelt 


endif 
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# 


c 


c 


c 


c 


c 

c 

c 

c 

c 

c 


else 

if (dmelt.gt. 0. 0) then 

if ( rmelt (i+l,j) .ge. 1.0) 

& rmelt (i,j) = rmelt (i, j)+dmelt 

else 

if ( (rmelt(i-l, j) .le. 0.0) .or. (qaz .It. 0.0)) 
& rmelt (i,j) = rmelt ( i , j ) +dmelt 

endif 
end if 

if (rmelt (i, j ) .gt. 1. 0) rmelt(i,j) = 1.0 
if (rmelt (i,j) .It. 0.0) rmelt (i,j) = 0.0 
qij = (romelt(i,j) -rmelt (i,j) ) *pcmmass*hf 1/dtime 
qlat = qlat-qij*dz 

abm = (-exp(-cl*(rpo(j) -rs) )+exp(-cl*(rpo(j) -rn) ) ) 
qr = qaz*abm + qt 
qij = qij + qr 

if (i .le. (ni-j+2)) then 
if (i.eq.2) then 

odmelt = (rmelt ( i , j ) -rpmelt ( i, j ))*. 5*dy* 

* (densol/denliq-1) 
else 

odmelt = (rmelt ( i , j ) -rpmelt ( i , j )) *dy* 

* (densol/denliq-1) 
endif 

ro(j) = ro(j) + odmelt 
endif 

if (rmelt ( (ni-j+2) , j ) -0 . 0 . le. 0 .000001) then 
ro(j) = ni*dy 
ne(j) = ni-j+2 
endif 


odmelt = 0.0 
40 a(i) = -as 

b(i) = an+as+ae+aw+apo+dapo+hi j 
c(i) = -an 

d(i) = (apo+dtapo) *to ( i , j ) +qi j+ae*te+aw*tw+hi j *tci j 
delapo = apoo(i,j) - apon(i,j) 
if (delapo. ge. 0.0) then 

qmove = delapo* ( 0 . 0-to ( i , j )) *dz 
else 

qmove = delapo* (0 . 0-t ( i , j )) *dz 
endif 

qsens = qsens+ (apo+dapo) *dz*t (i , j )- (apo+dtapo) *dz*to ( i , j ) 
qrad = qrad + qr*dz 
if (i.ne.l) tac = tac+abm 
if (i .eq. ne(j)) goto 60 
50 i = i+1 
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goto l 
return 
end 
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subroutine property ( cv , den , i , j , k) 

common/ bloke/rme It (15,10) , denliq, densol , walls , cvmet , denmet , 
+ kmet , kliq, ksol , cvliq, cvsol 

real kliq, ksol , kmet , k, kpure 


kpure = 1.0 / ( rmelt ( i , j ) /kliq + (l-rmelt(i, j) )/ksol) 
cvpure = rmelt ( i , j ) *cvliq + ( 1-rmelt ( i , j )) *cvsol 
denpur = rmelt ( i , j ) *denliq + ( 1-rmelt ( i , j )) *densol 

(2‘k‘k'k-k'kic’k‘k-k'kic’k'k‘kic‘k‘kje'k‘k‘k-k-k’k'k‘k'klc'k'k'k‘k'k-k-k'kjc'k'k‘k-k'k'k **★★**★***★★**★**★** 

c put storage material in parallel with can material * 

C*************************************************************** 

k = kpure* ( 1-walls) +kmet*walls 
den = denpur* ( 1-walls) +denmet*walls 

cv = (cvpure*denpur* (1-walls) +walls*cvmet*denmet) /den 

return 

end 



subroutine solver (i) 

common/bloka/a ( 15 ) ,b(15) , c(15) ,d(15) 
common/blokb/ni,ne(10) , j ,t(15, 10) 
m = ne ( j ) 

do 1 i = l,m-l 

r = a(i+l)/b(i) 
b ( i+1) = b ( i+1) -r*c ( i) 
d ( i+1) = d ( i+1) -r*d ( i) 
continue 

t(m,j) = d(m)/b(m) 

do 2 i = m-1, 1, -1 

t(i,j) = (d(i)-c(i)*t(i+l, j) )/b(i) 
continue 


return 

end 



subroutine heatgas(nj) 
real ntu 

common/blokb/ni , ne ( 10) , j , t (15 , 10) 

common/blokc/ thmelt , tgasin , gasmas , cpgas , ntu , del tp , 

+ deltdz 

if (j .eq. 1) then 
qgasp = 0.0 
qgasdz = 0.0 
end if 

temp= (t ( 1 , j ) +thmelt) -tgasin 

qgasp=qgasp+gasmas*cpgas*temp* ( 1-exp ( -0 . 5*ntu) ) / (nj ) 
qgasdz=qgasdz+gasmas*cpgas*temp* ( 1-exp (-ntu) ) / (nj ) 

if (j .eq. n j ) then 

deltp=qgasp/ (gasmas*cpgas) 
deltdz=qgasdz/ (gasmas*cpgas) 
endif 

10 continue 
return 
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